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Frequently Used Nomenclature 


Note: the units used in this work are such that c = 1 unless otherwise noted. Also, 
as described in section 3, any function of energy can be written as a function of range 
instead. 


a j( E ) 


inverse of mean free path of particles of type j and energy E ; 


its dimensions are 


length 


a,jk(E,E') macroscopic spectral distribution for production of particles of type j and 

energy E from particles of type k and energy E'\ 
its dimensions are ; — -A 

length -energy 

a,jk(Q, O', E, E') macroscopic cross section for production of particles of type j moving- 
in direction 0 with energy E from particles of type k moving- 
in direction O' with energy E': its dimensions are 1 

microscopic spectral distributions 


da 

dE 


E 

'//!•''• r /) 


length-energy-solid angle 


kinetic energy; unless stated otherwise, all energies in this work are kinetic 

defined by equation (22) or equation (34); it represents a 

“source” of particles of type j and range rj 

range of particles of type j, defined by equation (19) 


range of protons 



Sj(E) 

Uj{E) 
u jk {E,E') 
Ujk{Q, fi', E, E ') 

°i(E) 

a jk {E,E') 
a jk (n,ft,E,E') 
, E) 

c fijC at, o, e) 


the stopping power of particles of type j and energy E- 
it has dimensions of 

length 

same as a t (E). except it only includes decays 

same as cij k (E,E') except it only includes production from decays 

same as a,j k (VL,VL' ,E,E') except it only includes production from decay 

same as cij(E ), except it only includes collisions 

same as cij k (E,E') except it only includes production from collisions 

same as 0^(0, Q', E, E') except it only includes production from collisions 

flux of particles of type j located at position ~Ht\ traveling with 

energy E\ its dimensions are 

flux of particles of type j located at position traveling with 


velocity it = vQ and energy E: its dimensions are 


l 


area-solid angle-energy 


IV 



1 Introduction 


To design a spacecraft with optimum shielding, the ability to obtain detailed knowledge 
of the radiation in the craft’s internal environment, when given the external radiation 
conditions and a set of shielding specifications, is required. The NASA Langley Research 
Center code HZETRN has been developed for this purpose. This code does not yet include 
the effects of pions or muons, which are thought to partially account for the discrepancy 
between HZETRN predictions and a recent experiment [1]. The intent of this paper 
is to present a method and its implementation as a set of fortran subprograms, called 
MESTRN, for the inclusion of pions and muons in the NASA transport codes. Results 
for transport of a cosmic ray spectrum through aluminum and water are also presented. 

The development of a transport code for muons and pions was modeled on the existing 
ion transport methods [2], A perturbative solution to the Boltzmann transport equation 
(using the straight-ahead approximation) was used. This solution differs from the solution 
used for light ions only by making no assumptions about the spectra of produced particles. 
A subroutine that iterates this solution numerically, so that transport over large distances 
can be achieved, has been partially integrated into HZETRN. Using this modified code, 
some results have been obtained for transport through aluminum and water. 

In addition to the inputs normally required to run HZETRN, cross sections and stop- 
ping powers that describe muon and pion reactions are also required. Muon and pion 
stopping powers were scaled from proton stopping powers by using a scaling factor deter- 



mined by the Bethe theory of stopping power [2] . The production of pions was assumed 
to be exclusively from nucleus-nucleus and pion-nucleus reactions. A combination of for- 
mulas from Ranft [3] and Blattnig et al. [4] was used for parameterizations that describe 
cross sections for these reactions. Pion decay was assumed to be the only method of muon 
production. Because this decay is a two-body decay, the rate of muon production can be 
obtained by conservation of momentum and a Lorentz transformation. 

The content of this paper is as follows: Section 2 reviews the stationary Boltzmann 
transport equation and the continuous slowing down and straight-ahead approximations. 
Section 3 develops an analytic perturbative solution to the transport equation that is valid 
for transport over short distances. The latter part of this section shows how the analytic 
solution can be iterated numerically to achieve transport over arbitrary distances. Section 
4 contains the interaction database used as input into the transport equation, including 
the various pion and muon production cross sections. Section 5 presents results for charged 
pion and muon fluxes after the initial cosmic ray fluxes were transported through various 
depths of aluminum and water. Some conclusions are presented in this section as well. 
Appendix A contains many details of the calculations of the analytic solution to the 
transport equation that were left out of section 3 for the sake of clarity. Appendix B 
explains some of the notations that were used for cross sections in this paper. Appendix 
C contains the details of derivations of macroscopic muon cross sections that were listed 
in section 4. Appendix D describes the MESTRN Fortran subprograms, and appendix E 
contains these subprograms in their entirety. 
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2 Derivation of the Transport Equation 


Consider a region of space filled with matter described by the appropriate cross sections. 
Further, consider a spherical volume element of radius 8 centered about position it. 
Define <pj{lt , Q, E) to be the flux of particles of type j located at position x, traveling 
with velocity if = vQ and energy E. The relation between flux and the distribution 
function f is discussed in Reif [5], flux is related to f by a factor of velocity. The particle 
fluxes are assumed to be time independent, including only stationary states. Neglecting 
time dependence, the units of the flux 6Ax ,Q.,E) are pt — \ — . 

Now, consider the flow of particles through this volume element. The number of type 
j particles leaving a surface element 8‘ 2 dQ of the sphere in a direction Q. with energy 
E , is S 2 dQ(/)j(lt + 5Q.il, E). Similarly, S 2 dQ(f)j(lt — 5Q,£l,E) is the number entering 
the opposite part of the sphere. Figure 1 gives a sketch of the situation considered here. 
The number of particles entering the sphere differs from the number leaving by the gains 
and losses due to interactions with the material and the decay of unstable particles. Any 
interaction of the radiation with itself is assumed to be negligible. These considerations 
can be represented in a mathematical formula as [2, 6] 



S 2 dh[(f) j {lt + 8Q : h, E) - Ijilt - hO, O, E)\ = 
[ dE'dtt' n', E, E')<t> k {lt + m, E') 


-s 2 dh / di<f>j(it + m,n,E) aj {E) 


(i) 
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Figure 1: Spherical volume element of radius S with fluxes of particles at opposing posi- 
tions on labeled surface. These two fluxes are equal if no interaction occurs while particles 
move length S through sphere. 

The limits of the integrals over dE' and dfl' are constrained by conservation of momentum 
and energy, or equivalently, dE' and dW range over all values, and the cross sections are 
constrained by conservation of momentum and energy. 

Radiation is treated as a flux of classical particles; interference effects between the 
radiation and the bulk medium are therefore neglected. This type of interference is not 
expected to have much of an effect on the processes considered here because of the large 
energies of the particles in question. Quantum effects are, however, the main considera- 
tion when calculating cross sections for interactions between the radiation and individual 
nuclei. 

Equation (1) can be understood as a statement of conservation of flux. The change 
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in the amount of flux (the left-hand side) is equal to the amount destroyed (second term 
on right-hand side) subtracted from the amount created (first term on right-hand side); 
<Pk{~^ + fi, E)cijk{£l, ST, E, E') is the number of particles per volume, energy, and 
solid angle of energy E (and moving in the direction 11) that are created by particles 
of type k moving in the direction Q 1 with the energy E' . This statement defines what is 
meant by the macroscopic production cross section aj k (Q,Q',E,E'), which has the units 

of ] — 7T- — tt-t r . An example of this type of process is a proton (particle k) colliding 

with the medium and producing a pion (particle j). Particles can be created either 
from the decay of unstable particles or from interaction with the medium. For example, 
pions can be created from the collision of heavy ions with a medium or by the decay 
of kaons. Multiple reactions can produce these types of results, making the appropriate 
cross sections inclusive. 

a jk {Q, O', E, E') = a jk {h, O', E, E') + u jk {h , O', E, E') (2) 

where a lk represents production from collisions, and Uj k represents production from decay. 

By integrating over E' and O' and summing over k, incoming particles of all energies, 
angles, and types are considered, leaving the number of particles of type j produced per 
volume. To get the number of particles produced, we integrate over the path through the 
sphere f^ g ell and multiply by the surface element 5 2 clfl (see fig. 2); aj(E)<p k (h t + dQ, Q. E) 
is the number of particles of type j moving in the direction Q with energy E per volume 
that are lost. What is meant by lost is that they are no longer of type j, or they are 


5 



no longer moving in the same direction with the same energy. Furthermore, cij(E) is the 
inverse of the mean free path and has the dimensions length -1 . It is a value of unity 
divided by the distance a particle of type j and energy E will travel before having an 
interaction or decaying and becoming “lost”; cij(E) can be split into two components, one 
for decays and one for collisions. 

a j(E) = Uj(E) + (Jj{E) (3) 

The term Uj(E)(pj(E ) is the number of decays per unit volume per unit energy, and 
a.j(E)(pj(E ) is the number of collisions per unit volume per unit energy. The symbol 
( 7j(E ) is the total macroscopic interaction cross section (inelastic plus elastic) for particles 
of type j that interact with the medium. For the sake of simplicity, these cross sections are 
considered independent of position; the medium is assumed to be constant. See section 4 
for more details on cross sections and decay lengths, particularly the relationship between 
microscopic and macroscopic cross sections. 

The symbol S is the radius of the region under consideration, and it can be made 
arbitrarily small, allowing the fluxes to be expanded in a three-dimensional Taylor series 
around the vector ~T\ 

0(1^ ± Ifl, Q, E) = 0(T\ Q, E) ± IQ ■ v0(^, fh E ) + 0(d 2 ) (4) 
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Q 

<f>(lt + IQ, Q' , E ') 

Figure 2: At point it + IQ, particles moving in direction Q' create particles moving in 
direction Q. 

Equation (4) is a three-dimensional Taylor series expansion, in which only terms up to 
order 5 are kept. The operator Q ■ y is the spatial derivative in the direction Q. For 
example, if Q is pointed in the r direction, then Q ■ y = J:- Equation (5) is the integral 
of the fluxes given by equation (4). The terms of odd order in l integrate to 0, so there is 
no term of order S 2 in equation (5). 

The next step in the derivation is the substitution of the above two equations into 
equation (1) and the simplification of the result. A three-dimensional transport equation 
results [2, 6], 

0j{lt, Q, E) + 5Q ■ \70j{!t, Q, E) - {^{it, Q, E) - SQ ■ \/u ; ( r. Q, E)] + 0{6 2 ) 
clE' dQ' ^ a jk {Q, Q', E, E')[ 26(/> k (lt, Q', E')\ 

k 

- [2S<J)j(lt, Q, E)\cij(E) + 0(d 3 ) 
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which reduces to [2, 6] 


Q • yo,( A ft, E) = I dE'cin' Y ci jk {Q, ft', E, E' )<>,.{ A ft', A) 

^ k 

— E)cij(E) + 0(5) (6) 


The limit as 5 — >■ 0 can now be taken, making the terms of order 5 vanish. 


2.1 Continuous Slowing Down Approximation 

Cross sections Oj k and <jj describe all the relevant processes except those that involve 
decay. A simplification can be made if the atomic and nuclear components of these cross 
sections are treated separately. Examples of atomic cross sections are atomic or molecular 
transitions of the medium caused by the radiation. The nuclear cross sections describe 
collisions between the nuclei of the incoming radiation and the nuclei in the material 
under consideration. More precisely, what is meant by “atomic” is any process that can 
be approximately described by equation (9), and “nuclear” refers to all processes not 
accounted for elsewhere. The separation is made as follows: 


a jk (h,Q',E,E') 


a 


* k (n,tt,E,E') + of k (n,tt,E,E') 


(7) 


and 


= *?( E) + af(E ) ( 8 ) 

where superscripts n and at refer to the nuclear and atomic components, respectively. 
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No significant deflection occurs in a collision with an atomic electron when the fol- 
lowing conditions are met: the particle is much more massive than an electron (mj W 
melee iron ) • and the energy transferred is much smaller than the original energy ( E e n ) 
and is below particle production threshold. When these conditions are satisfied, the fol- 
lowing approximate form for af k can be assumed [2, 6]. 

W, E, E') (E)S 2 (h • - l)S jk 5(E + e n - E') (9) 

n 

where n labels the excitation level, e n is the energy of the excited state n, and af n (E) is the 
cross section for excitation to the nth state. The term S(E + e n — E') enforces conservation 
of energy, and the term 8j k preserves particle type. The angular delta function conserves 
direction, and is defined by the equation f f(Q')S 2 (Q ■ Q' — 1 )clW = /(fi), where £2 • Q 1 is 
the scalar product of two unit vectors. 

For the sake of clarity, only the terms involving atomic cross sections from equation (6) 
will now be considered, ignoring terms including nuclear components of the cross sections, 
decay lengths, and spatial derivatives. The terms of interest from equation (6) are 

f cm'dE'af^^.E.E 1 )^ ,h',E') - af^E)^ ,h,E) 

k ■’ 

Equation (9) then implies that [2, 6] 

^2 j cltt 1 clE'af k {n , E, E')(f) k (-£, fi', E') - afjE)^^, fi, E) 

k 

= ^VjniE + ,fi, E + e n ) - a" j , {E)o j { A <*>. A'! (10) 
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A series expansion of af n (E + e n ), keeping terms up to order e n , can now be made [2, 6]. 


^2 a fn( E + ^ E + e n ) - n“ j l (E)n j ( /•.') = ^af n {E)(j) j (^ ,C2,E) 

n n 

+ Y^\o%(E)^Ae)} + 0(4)-af(E)Mt,tlE) 

n 

« A [S . (E)0 .(^S1,£)] (11) 

where the stopping power is defined as Sj(E ) = '^2 n e n a f n {E). The cross section for 
an arbitrary excitation <jf{E) is simply the sum of the cross sections for a particular 
excitation af n {E). Thus [2, 6], 

of(E) = Y,of n (E) (12) 

n 

because the first and third terms on the right-hand side of equation (11) have canceled. 
The term af n {E) is a macroscopic cross section and has units of length -1 , not of area. 
Stopping power therefore has units of and can be interpreted as an average energy 
loss per unit distance traveled. Keeping terms only up to order e n is justified because 
E e n , but in some situations, such as nearly mono-energetic beams, higher order terms 
can become important [2], 

The continuous slowing down approximation is equivalent to the situation where parti- 
cles continuously lose energy as they travel through a medium, hence the name. The rate 
of energy loss is averaged over all different possible excitation levels n, and is given by the 
stopping power Sj(E). See [2, 7, 8, 9, 10, 11] for more details on these approximations, 
and on stopping power in general. 
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Substituting equation (11) into equation (6) and rearranging the terms results in a 
three-dimensional time independent transport equation in the continuous slowing down 
approximation [2, 6], 


n • v - f^s, (E) + «_,(£)] 40, n, E) 


dE' (& V aj 0 ft, E, E')40 , ft, E') 


(13) 


k 

where all of the cross sections are now purely nuclear. The superscript n was dropped to 
simplify the notation. Also recall that the derivative acts on both Sj(E) and <pj{E), 
as written on the right-hand side of equation (11). 


2.2 Straight-Ahead Approximation 

When particles are produced mainly in the direction of the primary flux, a further sim- 
plification can be made. This simplification, commonly referred to as the straight ahead 
approximation, reduces equation (13) from three spatial dimensions to one spatial dimen- 
sion. 

The straight-ahead approximation is valid when the particle production cross sections 
are highly peaked in the forward direction, which is when there is little change in direction 
due to collisions. This criterion is expressed in the following equation: 

a jk (Q, ft', E, E') ^ a jk {E, E')6 2 {Q. ■ Q' - 1) (14) 

This equation is often accurate in high energy collisions, in which conservation of mo- 
mentum keeps particles moving mainly in the forward direction. Because the direction of 
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radiation flow is not being greatly changed by collisions, only a single spatial dimension 
is needed. Thus, 

(f){~3^, Q, E) — >■ (f>(x, E) (15) 


and 

(16) 

resulting in [2, 6, 12, 13, 14] 


d d 
dx BE 


Sj( E ) + a j( E ) <f) j (x,E)= / dE'^2a jk {E,E')(f) k {x, 


Recall that cij(E) = <Jj(E) + Uj(E). Thus, in the above equation and 
dimensions. In the code HZETRN, the “length” x has units of and 


cross section a has the corresponding units The relation between 
and microscopic cross sections is 


E') (17) 

a have the same 
the macroscopic 
the macroscopic 


& P^microscopic v-^/ 

for each target species, where p = Numbe ;°; p^ 1 ” . . 

The straight-ahead approximation simplifies the transport equation by reducing it to 
a single spatial dimension. A three-dimensional problem can now be solved by super- 
position. Because of linearity, a superposition of the solutions to the single dimensional 
equation (17) is the solution to the full three-dimension equation (13). However, some 
assumptions must be made on the boundary conditions, which can also introduce some 
error as is described in Wilson et al. [2], 
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3 Solution Method 


Equation (17) is an integro-differential equation. The first step in our solution method 
is to transform equation (17) into a pure integral equation, integrating to get rid of 
the derivatives. Second, a series of variable transformations are made, simplifying the 
resulting equation. A perturbative solution accurate only over short distances is then 
obtained, but it can be iterated numerically to obtain results for longer distances. 


3.1 Analytic Solution 


First, define the residual range [2, 6, 14, 15] 


rAE) 


/■'' clE' 

o SAP) 


(19) 


which implies 


chj_ = 1 

dE Sj(E) 


( 20 ) 


Physically, r ; represents the average distance a particle of type j would travel if only 
atomic interactions (stopping power) were considered. 

Now, define a new function [2, 6, 14, 15], 


4’jfarj) = Sj(E)<pj(x, E) (21) 

Physically r ip :l . which has units of represents the amount of energy lost at location 

x by particles of type j and range r.j [2], 
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where A qj fa zMiEzIhljl and a = cij(rj). Derivations of equations (23) and (24) 

can be seen in appendix A. 

Equation (24) is an approximate algebraic solution to equation (17) which is valid for 
small h. However, the solution is for the function 'W not <p (recall eq. (21)). The flux <p can 
be obtained by simply dividing by the appropriate stopping power. For particle transport 


14 



that involves distances greater than a small value h, equation (24) can be iterated in steps 
of h in a manner that will be discussed subsequently. 


3.2 Numerical Solution to Transport Equation 


The first step in the numerical solution is to make all quantities a function of energy per 
amu (atomic mass unit) rather than energy. Equation (17) then becomes [2, 6, 14], 


- J-JE S ’ (£)+%(£)]<«*>£) = / dE'Y,«ME,E')<p„{x,E') (25) 


where A } is the mass of the jth particle in amu, which is approximately equal to the 
atomic number in the case of nucleons. Now define 


Sj(E) = t S,(E) 


fj(E) = 


dE' 


(26) 

(27) 


■Jo Sj(E') 

Equation (25) becomes exactly the same as equation (17), except with S rather than S. 
The rest of the analysis follows through exactly the same as in section 3.1, except with S 
rather than S in equation (24). 

Using the Bethe theory for stopping power implies [2, 7, 13] 


Sj(E) « ^-S k (E) 

Vk 

(28) 

jfj{E ) « Vk?k(E) 

(29) 


where zy = -f- and is charge. However, the Bethe theory, and consequently the previous 
equations, are only accurate for limited energy ranges [2, 7, 8, 13]. The simple relation 


15 



between ranges allows equation (24) to be written as a function of the range of one type 
of particle. By choosing this particle to be the proton (v p = 1), equation (24) becomes 
, u „ x ^ .exp(ah) - 1 


i'j{x + h, r p ) ~ exp (—ah) < ibj(x , r p + Vjh ) + Qj(x , r p + Vjh ) 


+A q.j exp (ah)(— — + — 
L a a 2 a‘ 


where 

a = aj(r p ) (31 

A _ fi j (br P + ^)-g J ^-/M>) (32 

J h 

and where all stopping powers (S) are now (S). In other words, the function -w :l is now 

t.\j(x. /• j = Sj(E)(t>j(x,E) (33 

and the new source term is 


Qj(x,r p ) = Sj(r p )^2 dr k a jk ( r p , r k ) tl' k ( x , r k ) 

k ■’ 

= ^(r p )^ / dr p a jk(r P ,r' p )S p (r p )Mx,r' p ) 


When both a <C 1 and ah <C 1, the quantity exp(o/i)(- — V) + ^2 in equation (30) can be 


sensitive to numerical round-off error. By expanding the exponential term in an infinite 


series, it can be shown that 


. . s(h 1 \ 1 h? r 2 1 2 

cp(ah) - + — = — 1 + -ah + — (ah) ex; 

V tx d / (X Z L o 1 ^ 


This relation is exact, and the term on the right-hand side of the equal sign is insensitive 


to round-off error when a < 1 and ah < C 1. 
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To simplify the numerical implementation, make the further approximation in eq. (30) 


Qj{x, r + vh ) Qj{x,r) (36) 

This term represents particles produced in collisions at position x with a range r + vh. 
Evaluating this term at range r neglects continuous slowing down for particles produced 
until the step after it is produced. For example, a pion is produced somewhere between 
x and x + h. Slowing down is not accounted for until x + h. Also, make another approx- 
imation for A q which is a first order approximation to the derivative of the term q. This 
approximation is similar to the previous one. Evaluate all terms at the same r rather 
than r and r + vh. These approximations have little effect on the results because h is 
small, and they greatly simplify the numerical implementation. 

Because everything is a function of one range, equation (30) can be put on a single 
grid of position and range, then iterated numerically in steps of h. Enter the initial %p : j 
into the right hand of equation (30), and ipj (at distances h away) will be given by the left- 
hand side. Numerical interpolation can be used to calculate values of r tp :l that do not lie 
directly on grid points. This process can be repeated to achieve transport over arbitrary 
distances. To calculate the fluxes for a three-dimensional geometry, simply calculate the 
single dimensional solution (eq. (30)) along rays leading to the points where knowledge 
of the radiation flux is desired, summing the results. 

Equation (30) gives the particle fluxes at position x + h as a function of the particle 
fluxes at position x, allowing this equation to be iterated numerically to achieve transport 
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over arbitrary distances. To numerically iterate, do the following: 

(1) Set up a discreet distance-range grid for the fluxes to be evaluated at each grid 
point. Since the solution depends only on position and proton range, the same grid can 
be used for all particle types. 

(2) Enter initial fluxes for all range values at x = 0. 

(3) Solve for fluxes at x = h by using equation (30). This step will need to use a 
numerical interpolation routine because the particle fluxes will be needed at values of 
range not lying directly on grid points. 

(4) Repeat step 3 until the desired depth is reached. 

3.3 Calculating the Source Term Numerically 

Because the source term (eq. (34)) appears as a function of range in the transport 
equation, and particle production cross sections are never calculated as functions of range 
(note that range depends on the material), there can be some difficulty in numerically 
evaluating equation (34). In HZETRN, there is both a proton range array and a proton 
energy array. The range in the /tli slot of the range array corresponds to the energy of 
the /tli slot in the energy array so that 

p £max 

q.j(x. /•;) = Sj{rj) y / dr'Spir'jcijkir^r'jfafar') 

fc J E m in 
/'Smax 

= Sjin)^ dE'a jk (Ei, E')<f> k (x, E') (37) 

£ ** ^min 
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where the integral is presently over the entire energy grid. Because il’k{x, E') is only 
stored in memory at discrete energy points, a numerical interpolation routine is needed to 
calculate the fluxes in the evaluation of the integral. The integral is calculated by using 
a Gaussian integration routine. 

3.4 Limits of Integration 

The values for E m \ u and E majX are different for pions than they are for muons. Pions are 
assumed to be created entirely from collisions of nuclei on nuclei and pions on nuclei. 
The pion production threshold for proton-proton collisions is approximately 290 MeV. 
Nucleus-Nucleus collisions are scaled from proton-proton collisions, so the threshold for 
nucleus-nucleus collisions is iust 290 Pion production threshold in pion-nucleus 

reactions is approximately 170 MeV or 1545^^ [16] when neglecting elastic scattering 
and charge exchange reactions. The threshold is defined as the minimum energy required 
to produce a particle. E m ; n is the minimum energy needed to produce a particle with a 
particular energy. E m i n can be approximated as the kinetic energy of the projectile or 290 
whichever is greater. For pion production, there is, in principle, no limit on /-.j nax ; 
however, a transport code will only transport particles of a finite energy range. For the 
case of space radiation, the initial particle fluxes fall off sharply at high energies, so some 
cut off energy is chosen and particles above that energy are not transported. Because the 
fluxes above this energy will be quite small, the effect of this cut off should be negligible 
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for the purposes of radiation shielding. E max for the pion is simply identical to this cut 
off energy. 

Muons are considered to be created entirely from pion decay, so E m in and E max are the 
minimum and maximum energies a pion can have and still decay into a muon of energy 
Ei. Since pion decay is a two-body decay, the muon energy is uniquely determined by the 
pion energy as described in appendix C. To calculate these limiting values, consider the 
Lorentz transformation muon energy from the pion rest frame to the lab frame 

= 7 n{E* + ruin - ftp* cos 9) - m,, (38) 

where E lt is muon energy in the lab frame, and 9 is the angle between the direction of the 
pion and the muon, 7 = E ^ n m is the Lorentz factor, and ft = 1 — ^ is a speed. E* is the 

muon energy in the pion rest frame, and p* is the corresponding momentum. Remember 
that all energies are kinetic unless otherwise noted. Pion energy will be minimum when 
the pion moves in the same direction as the muon, when 9 = 0. Pion energy will be 
maximum when the pion moves in the opposite direction as the muon, when 9 = tt. So, 

E p = r/A E * + m,, ± ft n p*) - m„ 

= 7t r ( E* + rrin ) ± p* ~ m n ( 39 ) 

which implies 

[E„ + % - 7 ,(E' + »>V)] 2 = (7 1 - 1)(P') 2 (40) 


20 



which has two solutions 


Ik 


( E * + m fl ){E u + m,,) ± p*p„ 


m. 


(41) 


The solutions corresponding to the maximum and minimum pion momenta are 


E- 


max/min 


( E * + m v )(E u + m^) ± p*p tx 


m-n 


m. 


(42) 


4 Database 

For the transport of pions and muons, all relevant interactions involving these particles 
need to be considered. The quantities describing these interactions when using the straight 
ahead and continuous slowing down approximations are as follows: 

(1) The inverse mean free paths ci(r ) 

(2) Stopping powers S(r) 

(3) Macroscopic particle production spectral distributions a,jk(r,r ) for all processes 
that produce significant numbers of pions and muons 

(4) The initial pion and muon spectrum, which for cosmic rays is 0 
These quantities will now be discussed in detail. 
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4.1 Inverse Mean Free Paths 

The inverse mean free path of a particle can be expressed as the sum of the total macro- 
scopic cross section (elastic plus inelastic) and the inverse decay length. 


%( r ) = ojM + u A r ) 


(43) 


Particle loss due to decay is described by the function Uj(r ), where u ) {r)(p ) {r) is the 
number of decays, per unit volume, of particles of type j and range r. To calculate this 
function, simply look up the lifetime of the particle in question [11], and perform a Lorentz 
transformation from the rest frame of the particle to the rest frame of the medium. Thus 
the inverse decay length is [12] 


Uj{r) = 


(44) 


k{r)crf/{r) 

where 7 = - e T)+ ? "j c ^he Lorentz factor, when Ej is the kinetic energy of particle j, 


and ft = ~ c = \J 1 — ft’Er) gi ves the velocity of the particle, and Tj is its lifetime. The term 
c is written explicitly here and is not set equal to 1 because cr is the quantity usually 
tabulated [11]. Elsewhere, c = 1 is assumed. 

The macroscopic cross section <Jj{r) represents all nuclear processes by which particles 
of type j either change species, energy, or direction; a j(r) therefore consists of both the 
elastic and the inelastic cross sections for all nuclear interactions with the medium. The 
cross section is macroscopic, and it will consist of the sum of all microscopic interactions 
weighted by the relevant densities. For a medium of N different types of nuclei, each with 
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a number density p n (units ^ ) , 

N 

°,(r) = Y. A.[<““ C M + <“'(>')] (45) 

n — 1 

where a®|f stlc (r) and a“ elastlc (r) are the microscopic cross sections for collisions of particles 
of type j with nuclei of type n. Note that the discreet nature of the medium on the micro- 
scopic level is, in effect, averaged over by the introduction of the continuous macroscopic 
density p. 

Because muons only interact by the electroweak force, their short-range interactions 
with nuclei can be neglected, and their long-range electromagnetic interaction can be 
included in stopping power; therefore, 

~ 0 (46) 

For pions, the nuclear reactions to consider are of the type tt + A — >■ x. However, 
elastic scattering can be ignored here because the energies of interest are high enough 
so that elastic scattering has little net effect. The net result of elastic scattering at 
high energies is a small change in direction and energy which are negligible in practical 
calculations. At low energies, energy loss due to nuclear collisions can become important 
and are included in a nuclear stopping power in HZETRN [2, 7]. For the analysis here, 
pions and muons have very low flux at low energies, so nuclear stopping should have little 
effect and will not be discussed further. Only inelastic pion nucleus reactions will be 
considered. Inelastic cross sections scale approximately with the size of the nucleus, so 
they will be approximately proportional to the cross section for pion-proton cross sections, 
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where the constant of proportionality is approximately equal to the ratios of the areas. 

2 

The cross-sectional area of a nucleus is area ~ A 3 where A is the atomic number of the 
target nuclei and the cross section is in units of barns. So, 

cr W A ~ A°- 73 a wp (47) 

where the more accurate value of 0.73 determined by Ranft [3] was used instead of two- 
tliirds. Parameterizations of microscopic pion-proton cross sections are given in the fol- 
lowing sections. 

4.1.1 Parameterizations for Reaction 71 Ap —? x 

Data for inelastic cross sections were scarce, so elastic cross section data from Groom et 
al. [17] were subtracted from the above parameterization for the total cross section to 
get more inelastic data. These “calculated inelastic data” were used in conjunction with 
experimental data [18, 19] to parameterize the inelastic cross section. 

The total cross section (elastic + inelastic) parameterization is split into two separate 
parts. One part is for low pion momentum, and the other is for high momentum and was 
taken from Groom et al. [17]. The data used were also taken from Groom et al. [17]. 


4.1.2 Parameterization for Reaction 7 r + p — >■ x 


For P lab < 54.89 (GeV/c): 


7T++P _ 

G t — 


a 1 


a 4 


1 + [{x - a 2 )/a 3 ] 2 1 + [{x - a 5 )/a 6 y 


&J 


atani - — 
0.5 + 


7 r 


a 10 (48) 
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where x = log(Pi a b). Also, o,\ = 4.422, ci 2 = —1.267, = 0.6411, <14 = 1.915, a-, = 0.3248, 

a 6 = 0.4602, ci 7 = 2.567, a 8 = 0.7135, a 9 = 0.3089, and ai 0 = 0.5633. 

For P lab > 54.89 (GeV/c): 

af +p = X s e + 1 1 s~ m - Y 2 s~ m ( 49 ) 

where s = mi + + 2 .()rn p \Jpl <lh + rrY is the invariant mass squared, and where 

X = 11.883, li = 28.59, Y 2 = 5.90, e = 0.093, = 0.358, and = 0.560 [17], A 

comparison of this parameterization with data [17] is given in figure 3. 



Figure 3: Data [17] plotted against the parameterization given by equations (48) and (49) 
for total cross section for the reaction 7r + + p — >■ x. 
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4.1.3 Parameterization for Reaction tv p — >■ x 


For P lab < 1.283 (GeV/c): 

tv- j-o 7 r „ _ / x — b 2 \ 21 , r „ _ / x — 65 s 2" 

er t +p = 61 exp — 0.5( — ) + 64 exp — O.o( — ) 

L 63 J L 06 J 

, f. / \ 2l 1 , [. / X ^10 \ 2"] ^ / r r\\ 

+b 7 l+(-) +69 1+ ( — - — -) (50) 

L Og J L On 

where x = logPi ab , 61 = 48.3, b 2 = —1.27, 63 = 0.182, 64 = 33.9, 65 = 0.903, b$ = 2.04, 
b 7 = 26.8, 6 8 = 0.0782, b 9 = 16.5, &i 0 = -0.333, and b n = 0.0744. 

For P Jab > 1.283 (GeV/c): 

af +p = AV + Yis~ m + Y 2 s~ m (51) 

where s = mi + m 2 p + 2.0m p y / pj ab + m 2 is the invariant mass squared, and where 
X = 11.883, li = 28.59, Y 2 = 5.90, e = 0.093, //, = 0.358, and //, = 0.560 [17]. A 
comparison of the parameterization for the reaction 7r _ + p — >■ x with data [17] is given 
in figure 4. 


4.1.4 Parameterizations for Reaction tt ± p — >■ x: Inelastic 


TT + +P __ 

; inelastic 


1 + exp[(c 2 - x)/c :i 


/ .r C5 , 

c 4 exp — O.o ( ) 

L v ^ ’ J 


where .r = log(Pi ab ), Ci = 20.1, c 2 = 0.0653, c 3 = 0.197, c 4 = 1.18, c 5 = 0.869, and 
c 6 = 0.705. 

The low energy parameterization for the inelastic cross section for negative pions on 
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Figure 4: Data [17] plotted against parameterization given by equations (50) and (51) for 
total cross section for reaction tt~ + p — >■ x. 


protons is the same form as that for the total cross section, but the parameters are 
different. 

For P ]ab < 5.02234 (GeV/c): 


_K + P 

inelastic 


cl i exp 


x — (I 2 \ 2l 


-o.M 


(I 4 exp 


-0.5(^— ^) 2 

dp, 


-\~dj 


1 + 0 2 


-1 


(i <) 


", / X — dip \ 

I- 1 dii > 


-1 


where x = logPi a b, ( U = 31.9, cl 2 = —1.27, cl 3 = 0.176, r/4 = 23.5, r/5 = 1.32, cle 
cl 7 = 12.5, d 8 = 0.0745, d 9 = 7.73, r/ K) = -0.33, and d u = 0.0518. 

For Pi ab > 5.02234 (GeV/c): 


(53) 

2.37, 


cP +p = 

inelastic 


AV + Y lS ~ m 


+ Y2S~ m 


(54) 
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where s = m^ + m| + 2.0m p y / p 1 2 ab + m 2 is the invariant mass squared, and where A" = 10.9, 
l'i = 20.2, Y 2 = 3.9, e = 0.093, ij\ = 0.358, and r\ 2 = 0.560. These parameterizations are 
shown in figures 5 and 6. 



Figure 5: Parameterization given by equation (52) for total inelastic cross section for 
7r + -proton collisions is plotted against experimental data [18] and “calculated data.” The 
latter was calculated by subtracting elastic cross section data [17] from total cross section 
parameterization given by equations (48) and (49). 
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Figure 6: Parameterization given by equations (53) and (54) for total inelastic cross 
section for 7r“-proton collisions is plotted against experimental data [19] and “calculated 
data.” The latter was calculated by subtracting elastic cross section data [17] from the 
total cross section parameterization given by equations (50) and (51). 

4.2 Stopping Power 

The Bethe theory shows that stopping power of other particles can be approximately 
scaled from the stopping power of protons [2, 7] 

Sj = ^S p (55) 

where Z is charge and A is mass in amu. The function S p is calculated in HZETRN and is 
an input to the pion transport subroutine. Because a positive pion has the same charge as 
a proton, it seems reasonable that the electromagnetic interactions could be simply scaled 
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from that of a proton to correct for the difference in mass. There might be some concern 
whether this is the case for negative pions and muons as well. The difference between 
stopping powers of protons and anti-protons was measured by Agnello et al. [20]. There 
seems to be an important difference at low energies up to the order of 100 IveV. Pions 
and muons are expected to have a similar behavior, and these energies are below where 
any significant pion or muon flux is expected. See [2, 7, 8, 9, 10, 11] for a more detailed 
description of stopping power. 

4.3 Production Cross Sections 

Cross sections a,j k contain particle production caused by collisions with the medium and 
the decay of unstable particles. For collisions, the macroscopic production cross section 
cijk{r,r') can be expressed in terms of microscopic spectral distributions and densities of 
target nuclei. In the straight ahead approximation, there is no angular dependence. Thus, 

r') = a jk {r , r') + u jk {r , r') 

N j 

= pn ( r ’ r ' ) + Ujk ( 56 ) 

n— 1 

where Ujk{i\r')(pk{r') is the number of particles per unit volume of type k with energy 
E' decaying into a particle of type j and energy E. The reaction being considered is of 
the type k — >■ j +x, where x represents any combination of particles and crj nk (r,r') is the 
microscopic cross section for a particle of type k with energy E'(r') colliding with a nuclei 
of type n and producing a particle of type j and energy E(r). There are many different 
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ways that this production can happen, so the appropriate cross sections are inclusive. In 
other words, the reaction is of the type A: + n — >■ j + x, where x represents any allowed 
combination of particles. The rest frame of the medium is assumed, so that the nuclei 
are stationary when neglecting thermal motion. Thermal motion can safely be neglected 
when the energies of the radiations are much greater than the thermal energy of the 
medium. 


4.3.1 Production From Decay 


The term u^ir. r') takes into account all possible decay channels. In other words, 


Ujk{r,r') = ^2u* jk (r,r') (57) 

i 

The index i indicates a particular reaction of the type a — >■ b + c + d, in which a, 6, c, and 
cl are all individual particles. 

As a first approximation, one can assume that muons are created entirely from the de- 
cay of charged pions. Charged pions decay by the reaction — >■ /A + /^ with a branching 

ratio of 99.99 percent. Because this is a two-bodv decay, the macroscopic particle pro- 
duction cross section can be calculated exactly by using only kinematic constraints, as 
shown in appendix C and in Gaisser [21]. 


^/i7 r ( E n i ) 


I -A - A 7^ 


(58) 


where p* = ^ J (E *) 2 — and E* = m "* + 2 ? " ;i — — and where rn„ & 0 is the mass of the 

neutrino. Note that the cross section u At7r (£’ At , E n ) has the above form only in the energy 
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region allowed by conservation of energy and momentum as specified in appendix C. 

Production of pions from the decay of other particles (kaons etc.) is presently ne- 
glected. 


4.3.2 Production From Collisions 


°jk(r,r') = 


/da \ try) 


k- 111 \j-i X 


(59) 


Presently, pions are assumed to be produced solely by the collisions, 


A + A 1 

-+ 

-)- x 

(60) 

7 r + + A 

-+ 

7T + + X 

(61) 

TT~ T A 

-+ 

7 r + + x 

(62) 

7T + + A 

-+ 

TT~ + X 

(63) 

7 T~ + A 

-+ 

7 r~ + x 

(64) 


where A represents a nucleus of atomic number A. Protons and neutrons are presently 
being treated equivalently with A set to unity. 

Parameterizations of inclusive pion production in proton-proton collisions are given in 
Blattnig et al. [4], All pion production cross sections are based on a combination of these 
parameterizations and parameterizations from Ranft [3]. These parameterizations assume 
energy in GeV while HZETRN has energy in MeV/amu. For the reaction A+A' -+ tt ± +x, 


, 0.69 , 1 / 0.69 / dcr\ 


( — \ — A u.by^/u.o» / — _ ) 

V (IE ) A+A'— Ytt+x V clE ) p+p— > tt+x 


(65) 
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For reactions of type x + A -+ 7 r + .x, the pion projectile is treated like a proton rescaled 


by the inelastic cross sections. 

( ( ^ a ^ 

' ( I /'. ' TT+A— >7T+.E 


^inelastic 
u ttA ^ 0.73 
.-.inelastic 
a pA 


da \ 

dE ) p+p— ytt+x 


0.62 A 0 - 77 


da \ 

dE) p+p— ’r'K+x 


( 66 ) 

(67) 


where parameterizations from Ranft [3] were used for the inelastic cross sections. Note 
that since the reaction p + p -A- n + x does not contain an elastic channel, the cross section 
for ^ + A ^ ^ + x will not include the elastic part. Since elastic reactions are also 
neglected as losses, elastic scattering is completely neglected in the present approximation. 
As discussed in section 4.1, elastic scattering has little net effect for the energies considered 
here. The following are parameterizations of charged pion production from proton-proton 
collisions from Blattnig et al. [4], The symbols E p and E n are the proton and pion kinetic 
energies in the lab frame. 

The positively charged pion spectral distribution for the range 0.3-2 GeV is represented 
by the following equations: 


F 2 

F l 

( ^ 

V dE ) lab 


= exp 


C+E% + C 3 E? 

( c 5 + ^ + r 7 /•:;/' + c 9 e$» + c u e^e^ + c 14 in /•;„) 


F 15 E^ lb — + C 17 E ^ 18 exp (Ci 9 \[Ek + C 2 o y/Ep) 


( 68 ) 


with constants C, given in table 1. 

The positively charged pion spectral distribution for the range 2-50 GeV is represented 
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Cl = 2.2 x 10 -s 

C's = -1-75 

Ci 5 = 2.5 x 10 s 

C 2 = -2.7 

Cg = -29.4 

Ci 6 = 0.25 

C 3 = 4.22 X 1CT 7 

Cio = 0.0938 

Ci 7 = 976 

C 4 = - 1-88 

C'n = -24.4 

Ci 8 = 2.3 

C 5 = 22.3 

Ci 2 = 0.0312 

Ci 9 = -46 

C 6 = 1-98 

C'i 3 = 0.0389 

C 20 = -0.989 

^ = -0.28 

Ci 4 = 1-78 



Table 1: Constants for Equation (68) 


by the following equations: 


Hi = 4.5 x 10- 11 

By = -35.3 

Bi 3 = 60322 

D-2 = -2.98 

D 8 = 0.0938 

B 14 = 1.18 

D 3 = 1.18 x 10“ 9 

Dg = -22.5 

Bib = -72.2 

£>4 = -2.55 

Bio = 0.0313 

Diq = 0.941 

B s = 22.3 

Bn = 2.5 x 10 6 

to 

II 

O 

He = -0.765 

B 12 = 0.25 



Table 2: Constants for Equation (69) 


F 2 

F 1 

f ^ 

V (IE ) lab 


d l e^ + d 3 e?* 

exp + ^ + D 7 E° 8 + D,Ey) 


— D\\E, 


D , 


F, 


DuE^ 14 exp (Di 5 y/EV + F>i G Ep 17 ^ 


(69) 


with constants Di given in table 2. 

The negatively charged pion spectral distribution for the range 0.3-2 GeV is repre- 
sented by the following equations: 


F 2 

Fi 


da \ 
(IE ) lab 


G t E^ + G 3 E° 4 


exp -4 — -j=- + G 7 E^ 8 + GgE^ 10 ^j 


E, 


Gi 


Ve, 

G\ 2 -fr + G\z exp(Gi4\/ / G 7r ) 
-c 2 


(70) 


with constants G, given in table 3. 

The negatively charged pion spectral distribution for the range 2-50 GeV is represented 
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G'i = 1.06 x 10“ 9 

C„ = -1.5 

G'n = 0.25 

Gi = -2.8 

Gr = -30.5 

Gi 2 = 2.5 x 10 s 

G ’ 3 = 3.7 x 10 “ 8 

G'g = 0.0938 

Gi 3 = 7.96 

Ga = -1-89 

G 9 = -24.6 

Gi4 = -49.5 

G 5 = 22.3 

Gio = 0.0313 



Table 3: Constants for Equation (70) 


by the following equations: 


Hi = 2.39 x 10“ 10 

H 7 = -31.3 

Hi 3 = 60322 

H 2 = -2.8 

Ilg = 0.0938 

Hu = 1.1 

H 3 = 1.14 X 10“ 8 

Ilg = -24.9 

fli5 = -65.9 

Ha = -2.3 

Rio = 0.0313 

H 16 = -9.39 

H 5 = 22.3 

Hu = 2.5 x 10 6 

1-In = -1.25 

H 6 = -2.23 

H 12 = 0.25 



Table 4: Constants for Equation (71) 


F 2 


F l 



lab 


= //,/ f"Wi // 3 /:" 4 

= exp (h 5 + -^+H 7 E«* + H 9 E»^ 

= H U E f# ^ + H U E^ exp (i /,5 + H U] E^) 

^2 


with constants i/j given in table 4. 


( 71 ) 


5 Results 

By using the solution to the transport equation given in section 3.2 and the database 
given in section 4, charged pion and muon fluxes were calculated after the initial cosmic 
ray fluxes were transported through various depths of aluminum and of water. Aluminum 
is commonly used in the construction of spacecraft, and water is a very large component 
of human tissue. The initial fluxes were given by the model for the 1977 solar minimum 
described in Wilson et al. [6]. The results of this calculation are given by figures 7-16. 
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Figure 7: Pion and muon fluxes after transport through 5 ^ of aluminum for 1977 solar 
minimum. 

Figure 7 shows pion and muon fluxes after transport through 5^ 0 f aluminum. The 
muon fluxes are too small to be seen. Figure 8 is the same as figure 7 except that it 
shows pion and muon fluxes as a percentage of the total calculated radiation flux. The 
percentages are large for higher energies. The peak in the percentages is at a higher 
energy than the peak in the actual fluxes. The peak in the percentage is at an energy 
at which the pion fluxes are relatively small. These large percentages at high energy are 
due to the falling off of the ion fluxes at high energy. Most of the radiation is at lower 
energies where the pion fluxes are relatively small. 

Figures 9 and 10, and figures 11 and 12 are the same as figures 7 and 8 except that 
transport through 10^|4r and 20^ of aluminum is shown. Figures 13 and 14 show pion 
and muon fluxes after transport through 20 of aluminum and S- 22 ! of water. Water is 

1 ° cm z cm z 
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Figure 8: Percentage of total radiation flux due to pions and muons after transport 
through 5 of aluminum for 1977 solar minimum. 

used to simulate tissue, and 5^ is the approximate depth of the blood forming organs 
in a person [2], Figure 15 shows the sum of the charged pion and muon fluxes at various 
depths in the aluminum and water. These fluxes steadily increase as more matter is 
traversed. The increase is due to an accumulation of pions as more nuclear collisions 
occur. Muons are negligible for these small amounts of material. Figure 16 is the same 
as figure 15 except that the percentage of the total radiation flux is shown. 

Muon fluxes are negligible compared to the pion fluxes, and the fluxes were also 
graphed separately in figure 17 so that the shape of the distributions could be seen. The 
fluxes are so low simply because there is not enough time for a significant amount of pions 
to decay. Muons might be much more important for thicker shielding such as the martian 
atmosphere. 
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Figure 9: Pion and muon fluxes after transport through 10 of aluminum for 1977 
solar minimum. 

In summary, a one-dimensional deterministic muon and pion transport code was de- 
veloped. This code is capable of transporting pions and muons over arbitrary distances 
for a wide variety of materials when it is coupled to the ion and neutron transport code 
HZETRN. All the needed cross sections were either derived or taken from other sources. 
Example results for transport through aluminum and water were then presented and 
shown to be reasonable. The major shortcomings of this calculation are as follows: 

1. The nuclear fragmentation database does take the effects of pion production into 
account. Because pion production is treated independently from other fragments, the 
total energy from all fragments in a collision will not necessarily add up to the energy of 
the initial colliding particles, so energy will not be strictly conserved. 

2. The parameterizations of cross sections describing pion production were based 
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Figure 10: Percentage of total radiation flux due to pions and muons after transport 
through 10 of aluminum for 1977 solar minimum. 

mainly on high energy data. The physical mechanisms for pion production are different 
at low energy, so the cross sections based on high energy data will be somewhat inaccu- 
rate at low energy. Improvements in these cross sections are presently being developed 
by using a meson exchange model where pions are produced from baryon resonances. 
Further work improving this code is presently underway. 
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Figure 15: The sum of pion and muon fluxes after transport through various depths of 
aluminum and of water for 1977 solar minimum. 



Figure 16: Percentage of the total flux due to the sum of pion and muon fluxes, after 
transport through various depths of aluminum and of water for 1977 solar minimum. 
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Flux amu/(cm MeV year) 



Kinetic Energy MeV /a mu 


Figure 17: Muon fluxes after transport through 5 ^ and 20 ^ of aluminum for the 


1977 solar minimum. 
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A The Solution of the Transport Equation 


As shown in section 2 and in equation (22), the continuous slowing down and straight- 
ahead approximations can reduce the transport equation (6) to the following equation: 


-d_ 

-dx 


_d_ 

drj 


+ 'MO! 




(72) 


where 

q j {x,r j ) = S j {r j )^2 / / d'r k a jk {r j ,r k )^ k {x,r k ) (73) 

k ■’ 

represents the creation of particles of type j, at position x, with range rj. The method 
of inverting this equation, and approximating the solution that was mentioned in section 
3 will now be discussed in detail. A similar analysis can also be seen in Wilson et al. [2] 
and in Shinn et al. [14], 

Equation (22) or (72) cannot be directly integrated because it contains two partial 
derivatives. If the following variable changes are made (method of characteristics), then 
the partial derivatives will become an ordinary single derivative, which can be directly 
integrated [6, 15]. 


> 1 ., O 

0 = + 0 ( 74 ) 
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therefore, 


i) i) dr/j i) d i) dr/j d d£j d 

dx drj dx drjj dx d dr.j dr/j drj d£j 

d d d d 

d>U + <Ki + (>>li <Ki 

= 2 jr- (75) 

Substituting equations (74) and (75) into equation (22) or (72) results in 

+ (76) 

This equation can now be integrated by introducing the following integrating factor [2, 6, 

14], 

= ex P 


1 r i ' 

2 <W-H 

J -Zj 


(77) 


which implies 


d <) f Vj £ ■ — //' 

2 % /,/( " J ‘ // ' /) = ./ dr} ' a j(' J ~2~ ) 

= l l A\ r >lj)«.j C i .. ;// ) 


(78) 


therefore [15] 


‘ d Vj 




£// I 'b S/ '// 


)/0(s/- / //i = 


\i ■ >U Zi >li 


) 


(79) 


A simple integral equation can now be formed by using equation (79) in equation (76) 
and integrating over ry' from — to , 


_d_ 

j-ii H 




■V/ I '/ S, 


)»s(.Zj,Vj) d v'j 


f" , , i n'i Zj n' is , r ,, 

/ d VjQj ( — 2 — ’ — 2 — W 


(80) 
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Through further simplification and the identification /i ? ( . — £,■) = 1, the following result 
can be obtained: 


giving [15] 


Vh (— ^ 2 Vi) ~ Vh(o, £,)/'( -£j) 

; , ,, ,, ,£j + Vj Zj n'j N 

d-Vj^Mv >l;ihi ( 2 • — 2 — ) 

S/ 




+ / dTfj»j{£j,Vj)Qj f 3 0 

■' -ii 


Vj ) 


(81) 


2 2 

Changing variables back to a; and r by inverting equations (74), explicitly writing out /a 
and simplifying notation by using 


/ //x /■'' ' r i >lj, 

a j (Vj ) = % ( g ) 

QjiVj) = Qi{ 2 • 2 > 

yields the following after some further algebra [2, 14], 

1 pX-Ti 


(82) 


( •/(./'. /';) = exp 


dVidjiVi ) • i\j) 


l r x ~ r i 

-- / dr/j exp 
^ J — x— r 7 


<Wji‘AVj) QjWj) 


(83) 


Writing out equation (83) in full using the definitions given in equation (82) results in the 
following, which is equation (23): 


Vj) = exp 


„ ,X + rj — T)! I 

d Vj ftj ( 3 ) i'3 (0, + rj ) 


1 

2 


dih exp 


1 

2 


, „ r '' ! O >l'L] (X + rj + rjfj x + rj-rfi 

d vM 5 -) M 3 -• 3 ) (S I! 
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A.l Simplification Through Variable Transformations 


Equation (83) or (23) can be further simplified by making a substitution in 

X + Tj + Tj” 

y = - 

J 2 

which implies 

2 cly' (86) 

2 y' — x — Tj (87) 

x + Vj — y' (88) 


dri- 


ll, = 


x 


r i ~ n. 


the first term, 


(85) 


The new integration limits can be set by noting that when V] * rj, then y' = x , and 
when rj" = — x — rj , then y' = 0. Equation (84) now becomes 




exp 

1 

+ 2 


Jo 


dy' a j(x + — y') ipj( 0, x + rj) 


dxj . exp 


dr/jCij ( 


x 


r J ~ Vj 


) 


x 


<h ( 


X 


Vj x I rj - rj'j 


(89) 


Making another substitution, y 


x + Tj + r/'j 
2 ’ 


results in 


O! 


exp 


dy' a j 


[x 


r j - y') 4>j(o,x 


L Jo 


r x 

"s? 

1 

^ c 
+ 

s~s 

•^5 

1 

i-H 1 C 

1 

1 

/ dy exp 

'0 

z J 2y—x—Vj z J 


X( lj (y,x + r j — y) 


(90) 
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a • i • / x + rj+ri" 

Again substitute y = ^ — -, resulting in 


i\j(x. Vj ) = exp — / dy'cij(x + rj — y') ^(0 ,x + rj) 


Jo 


dy exp 


Jo 


dy'cij (x + rj - y') qj ( y,x + rj-y ) (91 ) 


A. 2 Perturbative Solution 

If only transport over short distances is initially considered, equation (91) can be greatly 
simplified by making a series of approximations. Arbitrarily large distances can then be 
obtained by iterating the result. To simplify the notation, define 


cijiy') = cij{x + h + rj - y') 

Qjiv) = Qj{y, ^ + h + rj - y) (92) 

Now let x — >■ x + h in equation (91), where h is a small distance. See [2, 14, 15] for a more 
detailed discussion of what is meant by small. We have 


ipj(x + h, rj) = exp 


r-X + h 


JO 


x+h 


dy'cij (y 1 ) ii;,((),.r I // I rj ) 

rx+h . 


dy exp 


Jo 


dy' aj{y') '/,(//) 


(93) 
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This equation can now be simplified in the following manner. 


i'j(x + h,rj) = exp 


x px+h 

dy' + 

0 J x 

x px+h . 

dy + dy) exp 

JO Jx 1 

x+h 


dy'^j o-jiy 1 ) 


= exp 


dy'cijiy') 


exp 


ipj(0,x + h + rj) 

px+h . 

dy' + / dy'jcijiy') 


Qj(y) 


x exp 


= exp 


x+h px 

dy + ay j exp 
x Jo 

/ dy'a j (y l ) q 3 {y ) 
Jy 

x+h 


0 

x+h 


dy 1 a j(y') ?;•;((). x + h + Xj) 


dy'aj(y') 


dy'a j (y') 


exp 


dy'aj ( y ') ipj (0, x + rj + h ) 


a x+h px . _ px _ 

< ly+ J dy Jexp - J dy'ajiy') qjiy) 


(94) 


Recognizing that equation (91) implies 


4’j{x,Xj + h) = exp — / dy'aj(y') ipj(0^x + Xj + h) 


dy exp 


Jo 


dy'ajiy') n(y) 


(95) 


allows the two boldface pieces to be added, resulting in the following equation: 


il’j(x + h,Xj) = exp 


x+h 


dy'ajiy') si’jfarj + h) 


ax+h 


= exp 


dy exp[— / dy'cijiy'^q.jiy) 
Jy 

x+h 

dy'ajiy') i'A x ^ r j + h ) 


ax+h 


dy exp 


r-X + h 


dy’ajiy') Qjiy) 


(96) 
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Expanding this equation by using the definitions given in equation (92) yields 

px+h 

ipj{x + h, r.j) = exp — / cly'aj(x + h + rj — y ') 'ipj{x , rj + h) 


x+h 


cly exp 


rx+h 


dy’dj (z + h + rj - y’) qj{y , x + h + rj - y) (97) 


Now in the first term, make the substitution t' = y' — x, and in the second term, make 
the substitution t = y — x. Then, 


J i/’j{x + h , rj) = exp — / dt'cij{rj + h — t ' ) i'j(x , rj + h) 


JO 


r h 

- px+h. 

/ dt. exp 

— / dy'cij(x + h + rj — y) 

Jo 

J t+x 


(98) 


Make another substitution in the second term t' = y' — x. 


i)'j{x + h. rj) = exp — / dt'cij(rj + h — t' ) i rj + h) 


Jo 


dt exp 


Jo 


dt'cij(rj + h — t') qj(t + x, rj + h — t ) (99) 


Since h can be arbitrarily small, a series of approximations can be made to simplify 
equation (99) by expanding the integrand in a Taylor series and keeping terms up to first 
order in h. This approximation is similar to the approximations made in equations (4) 
and (5), 

dt'cij(rj + h — t ') ~ (h — t)aj(rj) (100) 

Since h is small, and t ranges from 0 to h in the integral, t is also small. The source term 
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qj(x + t, i'j + h — t) can therefore be approximated for small t. 


, , . . , . dqAx + t, r.j + h — t) . 

qj (x + t, r j + h-t) qj (x, r j + h) H - ^ \ t =o t 


dt 


qj{x , i-j + h) + A qjt 


( 101 ) 


where A q.j is an approximation of dq A x + t ^i +h with units of given by the fol- 

lowing equation: 

Aq . = '//(•'••'•./ i /; ) <I.M /; -Q) (|02) 

Using equations (100) and (101) in equation (99) and defining a = a.j ( r ? ) yields 
ipj{x + //. /',) = exp( ali)Cj(x. Vj + h) 


dtexp[a(t — li) <jj(x, rj + h) + A qjt] 


J o 


= exp (—ah) x {ipj(x,rj + h) 


dtexp(cit)[qj(x , r,- + h) + A qjt]} 


J o 


Finally, 


itj ( x + h. Vj) = exp {—cih) < r,- + h) + qj(x , r,- + /i) 


exp(a/i) — 1 


+A q~ 


h 1 \ 1 


exp(aft) (“ — 72 ) 


a a 2 / a 2 


(103) 


which is equation (24). 
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B Cross-Section Definitions 


There is no standard notation for cross sections. Below are some commonly used defini- 
tions of cross sections. The last ones are the convention used throughout this work. 

The total cross section is defined as [22, 23] 


^total — ^elastic 


^inelastic 


^absorption 


(104) 


The scattering cross section is defined as [2, 23, 24] 


^scattering — ^elastic 


^inelastic 


(105) 


giving 


^total — ^scattering T ^absoption 


(106) 


The reaction cross section is defined as [23] 


^reaction — ^inelastic 


^absoption 


(107) 


giving 


^total — ^elastic 


^reaction 


(108) 


Fernow [25] writes 


Cfiotal — ^elastic 


^quasielastic 


^absoption 


(109) 


Rolnick [26] writes 


Chotal — ^elastic T ^"inelastic 


( 110 ) 
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where (Jmeiastic = ^reaction, in other words, where all interactions except elastic scattering- 
are included in cr inelastic . The convention of Rolnick is used throughout this paper. 


C Two-Body Decays 


Recall that all energies are kinetic energies and that c = 1. The reaction to be considered 
here is of the type n — >■ fi + v. A general two-bodv decay can be treated simply by the 
appropriate replacements of 7 r, //,, and v. The production of muons will be treated here. 
Neutrinos will be considered only as they pertain to muon production. 

In a two-body decay, conservation of energy and momentum is enough to completely 
determine the kinematics. In the rest frame of the original particle a, production of 
particle b is isotropic, and its energy is [11] 


E’J 


ml + ml — ml 


2 m, 


— m, 


v- 


= E* 


( 111 ) 


In any other reference frame, energy is a function of angle that can be determined by a 
simple Lorentz transformation. 


E' = ^f(E + in + ip cos 0) — m 


( 112 ) 


The rate of decay is independent of the medium, so the production cross section is 
dependent only on the difference of the two angles Consider the following 
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equation: 


VilSV-Ej*. E,)fa(l£, E,)AE^(Pn^d 3 ~i = X'X.J fEi 


(113) 


where u, 17r ( fL, E,, ,E n ), which has the dimensions of i — -r- — tt- r is the macroscopic 

cross section for the production of muons from the decay of pions; d >] N ll7l . which has units 
ene 1 rgy , is the number of muons with energies between £), and £), + dE ^ moving between 
directions fp and fp + r/Q /( produced from pions of energy E % in the region between it 

and it + dlt. <p n (lt , E n ), which has the dimensions of , is the flux of pions of 

energy E n at position it . The cross section is independent of the direction of the pion, 
so only fluxes independent of pion direction need be considered. Now, also consider that 


N^{E n ) = 


flTT 


< ^ 7r -^tt) 


(in) 


where N fi7r (E n ) = f d (i N /J7r (E 7r ) and Ip is a branching ratio; r n is the mean lifetime of a 
pion, implying that /P-cy^r# is the mean decay length of a pion. Together equations (113) 
and (114) imply that 


r p, 7 r 

r 




I d 2 tl fl dE fi . 

'Ufl'K (O^, fp, E n ) 


(115) 


Ip = 0.999877 [11], which will be set to unity for the rest of this appendix. 

The production cross section of the muon in the rest frame of the pion can now be 
deduced from the previous equation. 


u 


(Q, ip, E v ) 


W ~ E*) 

An (31 s yl f rl f 


(116) 
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where rf signifies the rest frame of the pion. This quantity can be seen to have the correct 
angular dependence (i.e. , none), and the correct energy dependence, and it integrates to 
the correct result. The cross section is infinite since ft 1 ^ = 0 by definition. The transport 
equation needs cross sections in the rest frame of a medium, so a Lorentz transformation 
needs to be performed. 


C.l Lorentz Transformation 


Consider how the following quantity transforms under Lorentz transformations: 


^ 7T ( ) 


ft 77 TtT Tt7 


(117) 


Since equation (117) holds true in any reference frame, it can easily be deduced how 
u 7 T (E 1 r) transforms under Lorentz boosts. 


Equation (115) then implies 


<(K) = WZT^E,) 

Pttitt 


tTtt 


(118) 


n'^a'.ElE^dE^ia' = I u,A^E„,E,)dE,dn 

P'k'I'k 


7T ITT 

y 7T 


x'P'n 


= pipy / 

Pit In J 


d(W,E') 1 


(119) 


where = — is the Jacobian of the transformation between (O', £7) and (0,-E^), 

0{ll PfJ, r 1 

as shown in Hagedorn [27]. Hagedorn calculates the Jacobian for the total energy not 
kinetic, but clE toU \ = dEkmetic so the Jacobian is the same. 
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The macroscopic cross section in a boosted frame is 

e t . K) = El), E,(n\ £*); 

Pp Pit in 

The cross section in the lab frame is therefore 




5(E — E*) p 1 
47T/3 7r 7 7r T 7r p 

( s )[ 7 7r (£ ,/ + m + fi n p' cos 9') — E* — m]p' 
A-Kf3^Ti T \/i 2 {E' + m + fi Ti p' cos 9') 2 — m 2 


where all variables without subscripts from here on refer to the muon. 

In the straight-ahead approximation, a macroscopic spectral distribution is needed. 
Simply integrate over angle 


-^7r) J dQUjj, 71 - (£^5 -E'/i? r) 

= 2 n [ r/(eos 0)11,,, (Cl. ) 


Now make the substitution 


5[y$(E + m + fi w p cos 9) — E* — m = 


1 , E + m E* + m , 

%finP C0& + f%P JirfiirP 


resulting in 


^/i7T ( ^ 1 E , E v ) 


^s-S ( cos 9 + 

V PttP PttP ' 

47r/3 7r 7 7r T 7r \J ^ 2 {E + m + finP cos 9) 2 — m 2 


where the primes have been dropped for convenience. 

u^E,) = 2* Ucase) . „ ^ ^7 


J Aitfi^T* \/il(E + m + fi w p cos 9) 2 - m 2 

1 1 

AtTTtt 2fi^ v \J(E* + m) 2 — m? 

1 1 

^-^77 N — ItT 7t7 
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where E* is defined in equation (111) and p* = \J ( E * + m ) 2 — m3. Note that the cross 
section ?/, /J7r ( E lt . E n ) has the above form only in the energy region between 
7 (E* + m — j3p*) — rn and 7 (E* + m + j3p*) — m. It is zero outside this region due to 
conservation of energy. 

The total cross section, which is equal to the inverse decay length for the case of a 
single decay channel, should be the following: 

U/m ( E-k ) = — - — (126) 

Pit 7 Jr ^~7r 

therefore check that 


I dE^ Uf J/K ( Ej j , E w ) 

' E m in PirPir T ir 


(127) 


where E m in and E max are 7 (E*+m — {3p*) — m and 7 (E* + m+Pp*) — m, respectively. These 
limits are calculated by using the formula E = p{E* + m — f3p* cos 9*). The minimum 
energy occurs when 9 = 180° and the maximum occurs when 9 = 0°. 


"7 (E*+m.+/3 7T p*)-rn 


dE„ 


27t rArP* 


I r y(E*-\-m—(3 7r p*)—m P'k^'k^'k 




(128) 
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D MESTRN Subprogram Descriptions 


Unless otherwise stated: 


“Length” has the units 

All energies are kinetic and in units of LLY 

0 amu 


All fluxes are actually flux multiplied by stopping power as shown in equation (21) and 
have the units — . 


Stopping powers have the units ,cm2 

ir ir o i a,mn-erm 


Inverse mean free paths have the units 


9 

cm 
gm * 


2 

Inverse decay lengths have the units ^2-. 


Microscopic total cross sections have the units cm 2 . 


2 

Macroscopic total cross sections have the units £iS -. 

L orn 


Microscopic spectral distributions have the units cn L' A mi 


Macroscopic spectral distributions have the units 


All target properties needed for pion transport are included in the following common 


statement (see HZETRN user instructions and the main program of HZETRN): 
COMMON/TARGET/NLAY,XLAY,DN,NAT,ATRG,ZTRG,DENSTRG 
where 


NLAY is the layer number. 


XLAY is the thickness of a given layer. 
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DN is the density of the target in 


_gm_ 

cm 3 * 

NAT is the number of elements in the target. 

ATRG are the atomic numbers of the elements in the target. 

ZTRG are the numbers of protons in the elements in the target. 

DENSTRG is the density of the target in number . 

J O gm 

Function AAMMCS(EPROJ,EPI,P,J): This function calculates the macroscopic 

spectral distribution g v a(E v ,Ea) for production of pions of type .J with energy EPI from 
ions of type P with energy EPRO.J interacting with the medium. It has the following real 
input: 

EPRO.J is the energy of the projectile nucleus. 

EPI is the energy of the pion being produced, 

and the following integer input: 

P labels the type of nucleus projectile. 

J labels the type of pion produced. 

The function AAMMCS calls the following functions: 

AAPIP returns the microscopic spectral distribution for positive pion production in 
nucleus-nucleus collisions. 
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AAPIM is the same as function AAPIP except that it is for negative pion production. 


A returns the atomic number of ions. It is an entry in the HZETRN function TS. 

Function AAPIM(EPROJ,EPI,APROJ,AT): This function returns the microscopic 
spectral distribution || for the reaction A + A — >■ w~ + x. 

It has the following real input: 

EPROJ is the energy of the projectile. 

EPI is the energy of the produced pion. 

APROJ is the atomic number of projectile. 

AT is the atomic number of target. 

The function AAPIM calls the following functions: 

PIMINUS function returns the microscopic spectral distribution for reaction proton + 
proton — >■ 7 t~ + x 

Function AAPIP(EPROJ, EPI, APROJ, AT): This function returns the microscopic- 
spectral distribution || for the reaction A + A — >■ 7 r + + x. It has the same real input as 
the function AAPIM. 

The function AAPIP calls the following function: 

PIPLUS returns the microscopic spectral distribution for the reaction proton + proton — >■ 

7T + + X. 
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Function ANUMM(J): This function returns the dimensionless factor that scales the 

stopping power of pions and muons of type .J from that of protons. 

It has the following integer input: 

J is the index for pion/muon type (J = 1,2,3 or 4 correspond to // ' , jT . 7 r + , and tt~ , 
respectively). 

Function INELASTICPIMA(ATARG,EPI): This function returns microscopic to- 
tal (elastic + inelastic) pion-nucleus interaction cross sections a 7T (E) for the reaction 
7T _ T A — y x. It is temporarily set to return the inelastic cross section as discussed in 
section 4. 

Real input: 

AT is the atomic number of the target nucleus. 

EPI is the energy of the pion. 

Function INELASTICPIPA(ATARG, EPI): This function returns microscopic to- 

tal (elastic + inelastic) pion-nucleus interaction cross sections a 7T (E) for the reaction 
7r^~ -)- A — y x. It is temporarily set to return the inelastic cross section as discussed in 
section 4. 

Real input: 

AT is the atomic number of the target nucleus. 
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EPI is the energy of the pion. 


Function MAMMCS(EPROJ,EPI,P,J): This function returns macroscopic spectral 
distribution a 7T7T (E 7T , E w ) for production of pions from collisions of pions with the medium. 

It has the following real input: 

EPRO.J is the energy of the projectile pion. 

EPI is energy of the produced pion, 

and the following integer input: 

P labels the type of pion projectile. 

J labels the type of pion produced. 

The function MAMMCS calls the following functions: 

PIPAPIP returns microscopic spectral distribution for the reaction tt + + A — >■ tt + + x. 

PIMAPIP is the same as PIPAPIP except it is for the reaction tt~ + A — >■ 7r + + x. 

PIPAPIM is the same as PIPAPIP except it is for the reaction 7r + + A — >■ tt~ + x. 

PIMAPIM is the same as PIPAPIP except it is for the reaction tt~ + A — >■ tt~ + x. 

Function MGAUSSMM(A,B,N, SOURCE, J, EPI, E,PSI,PSIMMO, II, IJ,IM): This 
function is the Gaussian integration function mgauss slightly modified to integrate the 
function source over projectile energy. 
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It has the following real input: 


A is the lower integration limit. 

B is the upper limit. 

SOURCE is the function to be integrated. 

E is the energy array. 

EPI is the energy of pion/muon produced. 

PSI is the ion flux. 

PSIMMO is the pion/muon flux. 

The following integer input is used: 

J labels pion/muon type being produced. 

II is the number of points in the energy array. 

IJ is the number of ion types. 

IM is the number of pion/muon types. 

N is number of subintervals for the integration. 

The function MGAUSSMM calls the following function: 
SOURCE is the function to be integrated. 
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Function MMLOSS(EPI,J): This function returns the inverse mean free path for 

pion and muons of type J and of energy EPI. 

Real input: 

EPI is the energy. 

Integer Input: 

J is the particle type: 1 is p + ; 2 is jT: 3 is 7 r + ; and 4 is . 

The function MMLOSS calls the following functions: 

PIDECAY calculates the pion/muon inverse decay length. 

PNST returns the microscopic total (elastic + inelastic) pion-nucleus interaction cross 
section. 

Subroutine MULIMITS(EMU, THRESHOLD, UPPERLIM): This subroutine cal- 
culates integration limits for muon production from pion decay. 

Real input: 

EMU is the energy of the muon created. 

THRESHOLD is the lower limit. 

UPPERLIM is the upper limit. 
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Function MUPRODPI(EPROJ,EMU): This function returns the macroscopic spec- 

tral distribution w /J7r ( E lt . E n ) for production of muons from the decay of like charged pions, 
that is, for the decay tt ^ — )• /A + v^. 

Real input: 

EPRO.J is the energy of the projectile pion that is decaying. 

EMU is the energy of the muon created. 

The function MUPRODPI calls the following function: 

PIDECAY returns the inverse decay length of the pion. 

Subroutine OUTPUT: This function is temporarily used to output fluxes. It will be 

removed when pion/muon transport is fully integrated into HZETRN. 

Subroutine PGTMM (II,IJ,E,R,ST,PSI,PSIMM,H,QMOLD): This subprogram 

propagates the pion and muon components of radiation a distance h. 

It has the following real input: 

E is the energy array. 

R is the proton range array. 

ST is the proton stopping power array. 

PSI are ion fluxes. 
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PSIMM are pion/muon fluxes. 


H is the distance transported for one call of the subroutine. 

QMOLD holds array qm that was calculated during the last iteration 
and uses the following integer input: 

II is the number of energy/range grid points. 

IJ is the number of ions being transported. 

The function PGTMM calls the following functions: 

ANUMM calculates stopping power scaling factors for mesons. 

MMLOSS calculates inverse of the mean free path of pions and muons. 

PHI interpolates flux arrays. It is in the original HZETRN. 

SOURCE calculates pion and muon sources at specific energies. 

MGAUSSMM integrates the function source over projectile energy. 

Function PIDECAY(EPI,J): This function returns inverse decay length Uj(E ) of 

pions or muons. 

Real input: 

EPI is energy. 


70 



Integer input: 


J is the particle type: 1 is p + ] 2 is ^ i ; 3 is 7 r + ; and 4 is 7r . 

Function PIMAPIM(EPROJ,EPI,AT): This function returns microscopic spectral 
distribution || for the reaction n~ + A — >■ n~ + x. 

Real input: 

EPROJ is the energy of projectile pion. 

EPI is the energy of produced pion. 

AT is the atomic number of target. 

Function PIMAPIP (EPROJ, EPI, AT): Same as PIMAPIM but for the reaction 
tt~ T A — y 7r^ T x. 

Function PIMINUS(EPR,E): This function returns microscopic spectral distribu- 

tion for negative pion production from proton-proton collisions, that is, for the reaction 
p + p — >■ tt~ + x. 

Real Input: 

EPR is the energy of the projectile proton. 

E is the energy of pion. 
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Function PIPAPIM(EPROJ,EPI,AT): Same as PIMAPIM but for the reaction 

7T + + A — >■ TT~ + X. 

Function PIPAPIP(EPROJ,EPI,AT): Same as PIMAPIM but for the reaction 7r + + 

A — >■ 7T + + X. 

Function PIPLUS(EPR,EPI): This function returns the microscopic spectral dis- 
tribution || for positive pion production from proton-proton collisions, that is, for the 
reaction p -f p — >■ tt + + x. 

Real input: 

EPR is the energy of the projectile proton. 

EPI is the energy of the pion. 

Function SOURCE(EPI,EPROJ,E,PSIO,PSIMMO,J,II,IJ,IM): This function cal- 
culates pion and muon sources at specific energies. 

Real input: 

E is the energy grid. 

PSIO is the ion flux array. 

PSIMMO is the pion/muon flux array. 

EPI is the energy of the pion/muon produced. 
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EPRO.J is the energy of the projectile. 


Integer input: 

II is the number of range grid points. 

IJ is the total number of ions. 

IM is the total number of pion/muon types. 

J labels the type of pion/muon being produced. 

The function SOURCE calls the following functions: 

AAMMCS returns the macroscopic cross sections for production of pions from collisions 
of nuclei with the medium. 

MAMMCS returns the macroscopic cross section for production of pions from collisions 
of pions with the medium. 

MUPRODPI returns the macroscopic cross section for the production of muons from the 
decay of like charged pions. 

PHI interpolates fluxes. It is in the original HZETRN. 
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D.l 


Summary of Functions and Subroutines 


1. AAMMCS: macroscopic spectral distribution g^a^E^.Ea)- 

2. AAPIM: microscopic spectral distribution for A + A — >■ 7r _ + x. 

3. AAPIP: microscopic spectral distribution || for A + A — >■ 7r + + a:. 

4. ANUMM: factor that scales pion and muon stopping powers from proton stopping- 
powers. 

5. INELASTICPIMA: microscopic total cross section rr n [E) for 7r - + A — >■ x. 

6. INELASTICPIPA: microscopic total cross section a n (E) for tt + + A — >■ x. 

7. MAMMCS: macroscopic spectral distribution a 7T7T (E 7r , E v ). 

8. MGAUSSMM: integrates function SOURCE over energy. 

9. MMLOSS: inverse mean free path c ij(E ) of pions and muons. 

10. MUPRODPI: macroscopic spectral distribution u ll7T for tt ^ — >■ /A + v^. 

11. OUTPUT: temporary subroutine used to output fluxes. 

12. PGTMM: main propagation subroutine for pions/muons. 

13. PIDECAY: inverse mean decay length for pion/muon Uj(E). 

14. PIMAPIM: microscopic spectral distribution fj- for tt~ + A — >■ tt~ + x. 
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15. PIMAPIP: microscopic spectral distribution || for n~ + A — >■ n + + x. 

16. PIMINUS: microscopic spectral distribution || for p + p — >■ n~ + x. 

17. PIPAPIM: microscopic spectral distribution || for 7r + + A —> n~ + x. 

18. PIPAPIP: microscopic spectral distribution ^ for 7r + + A — >■ 7r + + :r. 

19. PIPLUS: microscopic spectral distribution || for p + p — >■ 7r + + :r. 

20. SOURCE: calculates pion and muon sources. 
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E The Meson and Muon Transport Code MESTRN 


THE SOFTWARE CODE PROVIDED HEREIN IS PROVIDED “AS IS” WITHOUT 
ANY WARRANTY OF ANY KIND, EITHER EXPRESSED, IMPLIED, OR STATU- 
TORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY THAT THE SOFT- 
WARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF 
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND FREEDOM 
FROM INFRINGEMENT, AND ANY WARRANTY THAT THE DOCUMENTATION 
WILL CONFORM TO THE SOFTWARE, OR ANY WARRANTY THAT THE SOFT- 
WARE WILL BE ERROR FREE. IN NO EVENT SHALL NASA BE LIABLE FOR ANY 
DAMAGES, INCLUDING, BUT NOT LIMITED TO, DIRECT, INDIRECT, SPECIAL 
OR CONSEQUENTIAL DAMAGES, ARISING OUT OF, RESULTING FROM, OR IN 
ANY WAY CONNECTED WITH THIS SOFTWARE, WHETHER OR NOT BASED 
UPON WARRANTY, CONTRACT, TORT, OR OTHERWISE, WHETHER OR NOT 
INJURY WAS SUSTAINED BY PERSONS OR PROPERTY OR OTHERWISE, AND 
WHETHER OR NOT LOSS WAS SUSTAINED FROM, OR AROSE OUT OF THE 
RESULTS OF, OR USE OF, THE SOFTWARE OR SERVICES PROVIDED HERE- 
UNDER. 

USERS OF THE SOFTWARE CODE PROVIDED HEREIN AGREE TO WAIVE 
ANY AND ALL CLAIMS AGAINST THE U.S. GOVERNMENT, THE U.S. GOV- 
ERNMENT’S CONTRACTORS AND SUBCONTRACTORS, AND SHALL INDEM- 
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NIFY AND HOLD HARMLESS THE U.S. GOVERNMENT AND THE U.S. GOVERN- 


MENT’S CONTRACTORS AND SUBCONTRACTORS FOR ANY DAMAGE THAT 
USER MAY INCUR FROM USER’S PRIOR OR FUTURE USE OF THE PROVIDED 
SOFTWARE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON, OR 
RESULTING FROM, THE USE THEREOF. 


The complete Fortran subprograms that make up the transport code MESTRN are 
contained in this appendix. In this code all energies are kinetic. Also, cross-section 
functions were named after the reactions. For example, AAMMCS returns cross sections 
for the reactions involving AA — >■ Meson/Muon + x. 

REAL FUNCTION AAMMCS (EPROJ , EPI, P, J) 

* This function calculates the the macroscopic spectral distribution 

* for production of pions/muons of type J from ions of type P. 

* 

IMPLICIT NONE 

* 
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* Integer Variables 

* NAT number of different types of atoms in medium 

* P labels type of ion projectile 

* J labels type of meson being produced 

* K loop index 

INTEGER NAT, P, K, J 

* Real Variables 

* DENSTRG density of atoms of type K in medium 

* ATRG atomic number of those atoms 

* EPROJ energy of projectile 

* EPI energy of meson being produced 

* AAPIP function that returns microscopic spectral distribution 

* for positive pion production in nucleus-nucleus 

* collisions 

* AAPIM same as AAPIP but for negative pion production 

* A function that returns atomic number it is an entry in 

* function TS 

REAL DENSTRG (5) ,EPR0J , EPI , SIGMA , A , AAPIP , AAPIM, ATRG (5) 


78 





* unused reals for common statement 

REAL NLAY , XLAY , DN , ZTRG (5) 

* set commons 

COMMON/TARGET/NLAY , XLAY , DN , NAT , ATRG , ZTRG , DENSTRG 

* 

* start function aammcs 

* 

* initialize aammcs 

* 

AAMMCS=0 . 0 

* 

* First consider positive pions 
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IF (J .EQ. 3) THEM 


* 

* 

* Slim micro-cross section over target atoms to get macroscopic cross 

* section 

* 

DO 10 K=1 ,NAT 

* 

* calculation based on scaling from pp collisions to AA collisions 

* 

AAMMCS= AAMMCS + DENSTRG(K)* 

$ AAPIP (EPROJ , EPI , A (P) , ATRG (K) ) 

* 

* end sum over target atoms 

* 
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10 


CONTINUE 


* 

* For negative pions 

* 

ELSE IF (J .EQ. 4) THEN 

* 

* sum micro over target atoms to get macroscopic cross section 

* 

DO 20 K=1,NAT 

* 

* calculation based on scaling from pp collisions to AA collisions 

* 

AAMMCS= AAMMCS+DENSTRG (K) * 

$ AAPIM(EPROJ,EPI,A(P) ,ATRG(K)) 

* 
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* end sum for negative pions 

* 

20 CONTINUE 

* 

* error statement for debugging 

* 

ELSE 

PRINT*, "error in aammcs there should be no else" 

* 

* end function AAMMCS 

END IF 
RETURN 
END 
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REAL FUNCTION AAPIM(EPROJ ,EPI , APROJ , AT) 


* this function that returns microscopic spectral distribution for 

* the reaction A + A -> pi- + x. 

* 

IMPLICIT NONE 

* 

INTEGER A 

* Real Variables 

* EPROJ kinetic energy of projectile nucleus 

* EPI kinetic energy of produced pi+ 

* PIMINUS function returns spectral distribution for reaction 

* proton + proton -> pi- +x 

* ALPHA scales collisions from target nucleon to nucleus 

* APROJ atomic number of projectile 

* AT atomic number of target 

REAL EPROJ, EPI, PIMINUS, ALPHA, APROJ, AT 
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* 

* start function AAPIM 

* 

* set alpha. 

* 

ALPHA =0.69 

* 

* calculate cross section 

* 

AAPIM=PIMINUS (EPROJ ,EPI) * (AT*APR0 J) **ALPHA 

* 

* end function aapim 

RETURN 

END 
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REAL FUNCTION AAPIP (EPROJ ,EPI , APROJ , AT) 

* this function returns microscopic spectral distribution for 

* the reaction A + A -> pi+ + x. 

* 

IMPLICIT NONE 

* 

* Real Variables 

* EPROJ kinetic energy of projectile nucleus 

* EPI kinetic energy of produced pi+ 

* PIPLUS function returns spectral distribution for reaction 

* proton + proton -> pi+ +x 

* ALPHA scales collisions from target nucleon to nucleus 

* APROJ atomic number of projectile 

* AT atomic number of target 
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REAL EPROJ, EPI, PIPLUS, ALPHA, AT, APROJ 


start function AAPIP 

* set alpha. 

* 

ALPHA =0.69 

* 

* calculate cross section 

* 

AAPIP=PIPLUS (EPROJ , EPI) * (AT*APR0J) **ALPHA 

* 

* end function aapip 
RETURN 
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END 


REAL FUNCTION ANUMM(J) 

* This function inputs the index for meson type (J = 1,2,3 or 4, 

* mu+, mu-, pi+, and pi- respectively) and returns the factor ANUMM 

* which scales the stopping power of mesons and muons to that of 

* protons. 

* 

IMPLICIT NONE 

* 

* Integer variables 

* J labels meson type 

* IM number of meson types 

* ZPI labels meson charge 

INTEGER ZPI, J, IM 
PARAMETER (IM=4) 
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DIMENSION ZPI(IM) 


DATA ZPI/1, -1, 1, -1/ 

* Real Variables 

* API is mass in amu 

REAL API (IM) 

DATA API /O. 11342, 0.11342, 0.14983, 0.14983/ 

* calculate 

* 

ANUMM= (ZPI ( J) ) **2/API ( J) 

* 

* end function anumm 

RETURN 

END 





REAL FUNCTION INELASTICPIMA (ATARG , EPI) 


IMPLICIT NONE 

REAL EPI, EPIGEV, PLAB, MPIGEV, API, X, A(ll), ATARG, 
$ GAUSSA , GAUSSB , LORA , LORB , S , B (6) ,MP 

DATA A /31.9, -1.27, 0.176, 23.5, 1.32, 2.37, 


$ 12.5, 0.0745, 

7.73, 

-0.33, 

0.0518/ 

DATA B /10.9, 20.2, 

3.90, 

0.093, 

0.358, 0.560/ 


MPIGEV=. 139570 
API=0. 14983 
MP=0. 939272 

* change units of energy to mev and calculate momentum 

EPIGEV=EPI*API/ 1000 . 0 

PLAB=SQRT (EPIGEV**2 . 0-EPIGEV*MPIGEV+2 . 0*MPIGEV**2 . 0) 

* calculate cross section 
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X=LOG(PLAB) ; 


IF (plab < 5.02234) THEM 

gauss a= A ( 1 ) * exp(-0.5*( (x-A(2))/A(3))**2.0) 
gaussb= A(4) * exp(-0 . 5* ( (x-A(5) ) /A (6) ) **2 . 0) 
lora=A(7) / ( 1 . 0+ (x/A (8) ) **2 . 0) 
lorb=A(9)/(1.0+((x-A(10))/A(ll))**2.0) 
IMELASTICPIMA=gaussa+ gaussb + lora +lorb 
ELSE 

s=mpiGEV*mpiGEV + mp*mp +2 . 0*sqrt (plab*plab 
$ + mpiGEV*mpiGEV) *mp; 

INELASTICPIMA=B (1) *s**B (4) + B (2) *s** (-B (5) ) 

$ + B(3)*s**(-B(6)) 

END IF 


* change cross section units from mb to cnT2 and scale to 

* pion nucleus collision 


INELASTICPIMA=INELASTICPIMA*1E-27*ATARG**0 . 73 
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* end function inelasticpima 

RETURN 

END 

REAL FUNCTION INELASTICPIPA (ATARG , EPI) 

IMPLICIT NONE 

REAL EPI, EPIGEV, PLAB, MPIGEV, API, X, A (6), ATARG 

DATA A /20.1, -0.0649, 0.197, 1.18, 0.866, 0.712/ 
MPIGEV=. 139570 
API=0. 14983 

* change units of energy to GeV and calculate momentum 

EPIGEV=EPI*API/ 1000 . 0 


91 



PLAB=SQRT (EPIGEV**2 . 0-EPIGEV*MPIGEV+2 . 0*MPIGEV**2 . 0) 


* calculate cross section 

X=L0G(PLAB) ; 

INELASTICPIPA=A(1)/(1.0 + exp( (A (2) -x) /A (3) ) ) 

$ + A (4) *exp(-0 . 5* ( ( (x-A(5) )/A(6))**2.0)) 

* change cross section units from mb to cnT2 and scale to 

* pion nucleus collision 


INELASTICPIPA=INELASTICPIPA*1E-27*ATARG**0 . 73 


* end function inelasticpipa 

RETURN 

END 
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REAL FUNCTION MAMMCS (EPRO J , EPI, P, J) 

* this function returns macroscopic cross section for production 

* of mesons from collisions of mesons with medium. Production from 

* decay is not included. 

* 

IMPLICIT NONE 

* 

* Integer Variables 

* NAT number of different types of atoms in medium 

* P labels type of meson projectile 

* J labels type of meson being produced 

* K loop index 

INTEGER NAT, P, K, J 
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* Real Variables 


* ATRG atomic number of those atoms 

* DENSTRG density of atoms of type K in medium 

* EPROJ energy of projectile 

* EPI energy of meson being produced 

* PIPAPIP function that returns microscopic spectral distribution 

* for the reaction pi+ + A -> pi+ + x 

* PIMAPIP same as PIPAPIP but for the reaction pi- + A -> pi+ + x 

* PIPAPIM same as PIPAPIP but for the reaction pi+ + A -> pi- + x 

* PIMAPIM same as PIPAPIP but for the reaction pi- + A -> pi- + x 

REAL DENSTRG (5) , EPROJ , EPI , PIPAPIP , PIMAPIP , 

$ ATRG(5) , PIPAPIM, PIMAPIM 

* unused reals for common statement 

REAL NLAY,XLAY,DN,ZTRG(5) 

* set common variables 

COMMON/TARGET/NLAY , XLAY , DN , NAT , ATRG , ZTRG , DENSTRG 
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* 

* start function 

* 

* initialize mammcs 

* 

MAMMCS=0 . 0 

* 

* First consider production of positive pions 

* 

IF (J .EQ. 3) THEN 

* 

* sum micro-cross section over target atoms to get macroscopic 

* cross sections 

* 
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DO 30 K=1,NAT 


* 

* calculate microscopic cross sections. 

* kaons are not included so no source from decay. 

* for the reaction pi+ + A -> pi+ x 

* 

IF (P .EQ. 3) THEN 

MAMMCS= MAMMCS + DENSTRG(K)* 

$ PIPAPIP(EPROJ,EPI,ATRG(K)) 

* 

* for the reaction pi- + A -> pi+ x 

* 

ELSE IF (P .EQ. 4) THEN 

MAMMCS= MAMMCS + DENSTRG(K)* 

$ PIMAPIP (EPROJ ,EPI , ATRG(K) ) 

* 
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* end sum over target atoms 

* 

END IF 

30 CONTINUE 

* 

* Now consider production of negative pions 

* 

ELSE IF (J .EQ. 4) THEN 

* 

* sum micro-cross section over target atoms. 

* 

DO 40 K=1,NAT 

* 

* calculate microscopic cross sections 
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* Also kaons are not included so no source from decay 

* for the reaction pi+ + A -> pi- x 

* 

IF (P .EQ. 3) THEM 

MAMMCS= MAMMCS + DENSTRG(K)* 

$ PIPAPIM (EPRO J , EPI, ATRG(K) ) 

* 

* for the reaction pi- + A -> pi- x 

* 

ELSE IF (P .EQ. 4) THEM 

MAMMCS= MAMMCS + DEMSTRG(K) * 

$ PIMAPIM(EPROJ ,EPI , ATRG(K) ) 

* 

* end sum 

* 
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END IF 


40 CONTINUE 

* 

* end function mammcs 

END IF 
RETURN 
END 

REAL FUNCTION MGAUSSMM ( A , B,N, SOURCE , J , EPI , E , FLUX , FLUXMM , II , 

$ IJ.IM) 

* This function is the Gaussian integration function mgauss 

* modified to integrate the function source over projectile energy 

* 

IMPLICIT NONE 

* 
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C function MGAUSSMM USES A GAUSS-LEGENDRE QUADRATURE FORMULA TO 
C PERFORM A 5 POINTS PER SUBINTERVAL NUMERICAL INTEGRATION 

* Integer Variables 

* N number of subintervals 

* XI =K-1 

* M,K loop index 

* J labels meson type being produced 

* II number of points in energy array 

* IJ number of ion types 

* IM number of meson/muon types 

INTEGER N, J, M, K, XI, II, IJ, IM 

* Real Variables 

* A lower integration limit 

* B upper limit 

* SOURCE function to be integrated 

* E energy array 

* EPI energy of meson/muon being produced 
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* EPROJ energy of projectile-variable of integration 

* FLUX ion flux 

* FLUXMM meson/muon flux 

* SUM result of integration 

C U(I) ARE THE ZEROS OF LEGENDRE POLYNOMIAL OF DEG. 5 

C G(I) ARE THE WEIGHTS SUCH THAT THE INTEGRATION IS EXACT 

REAL A , B , SOURCE ,EPI ,E (II) 3 FLUX(II 3 IJ) , FLUXMM (II 3 IM) 3 EPR0J 3 
$ FOFX ,U(5) ,G(5) 3 FINE 3 DELTA 3 UU 3 SUM 
EXTERNAL SOURCE 

* 

* start function 

* 

SUM=0 . 

IF(A.EQ.B) THEN 
MGAUSSMM=0 
RETURN 
END IF 
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C U(I) ARE THE ZEROS OF LEGENDRE POLYNOMIAL OF DEG. 5 

U(l)=. 42556283050918 
U(2)=. 28330230293537 
U(3)=. 160295215850488 
U(4)=. 067468316655508 
U(5)=. 013046735741414 

C G(I) ARE THE WEIGHTS SUCH THAT THE INTEGRATION IS EXACT 

G(l)=. 147762112357376 
G(2)=. 13463335965499 
G(3)=. 109543181257991 
G(4)=. 074725674575290 
G(5)=. 033335672154344 
FINE=N 

DELTA=FINE/(B-A) 

C FIND INITIAL VALUE FOR INTERVAL 
DO 40 K=1,N 
XI =K-1 

FINE=A+XI/DELTA 
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C EVALUATE FUNCTION AT 5 POINTS PER INTERVAL 

DO 20 M=1 , 5 

EPRO J=U (M) /DELTA+FINE 

FOFX=SOURCE(EPI, EPROJ, E, FLUX, FLUXMM, J, II, IJ, IM) 
C ADD TO SUM FOR EACH FUNCTION 

SUM=G (M) *FOFX+SUM 
20 CONTINUE 

DO 30 M=1 , 5 

EPRO J= ( 1 . -U (M) ) /DELTA+FINE 

FOFX=SOURCE (EPI , EPROJ, E, FLUX, FLUXMM, J, II, IJ, IM) 
SUM=G (M) *FOFX+SUM 
30 CONTINUE 

40 CONTINUE 

MGAUS SMM= SUM/DELT A 
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RETURN 


END 

* end function mgaussmm 

REAL FUNCTION MMLOSS(EPI, J) 

* This function calculates total interaction cross section + decay 

* length for mesons and muons of type J and of energy EPI 

* 

IMPLICIT NONE 

* 

* Integer Variables 

* J is the meson type: 1 is mu+ , 2 mu-, 3 pi+, 4 pi- 

* NAT number of different types of atoms in medium 

* K loop variable 

* 

INTEGER K, J, NAT 
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* 

* Real variables 

* ATRG atomic number of those atoms 

* EPI is meson energy 

* TOTAL is meson total interaction cross section with medium 

* PIDECAY is a function that calculates meson inverse decay length 

* DENSTRG density of atoms of type K in medium 

* PNST function returns microscopic total (elastic + inelastic) 

* pion-nucleus interaction cross section 

REAL EPI, TOTAL, PIDECAY, DENSTRG (5) , ATRG (5) , 

$ INELASTICPIPA , INELASTICPIMA 

* unused reals for common statement 

REAL NLAY , XLAY , DN , ZTRG (5) 

* set commons 

COMMON/TARGET/NLAY , XLAY , DN , NAT , ATRG , ZTRG , DENSTRG 
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* initialize total 

* 

T0TAL=0 . 0 

* 

* calculate macroscopic total (elastic + inelastic) interaction 

* cross section 

* First consider pions. charge differences between pions are ignored 

* 

IF (J .EQ. 3) THEM 

* 

* sum micro-cross section over target atoms to get 

* a macroscopic cross section 

* 

DO 10 K=1 ,NAT 

* 
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* calculation based on scaling from pi+ p collisions to pi+ A 

* 

T0TAL= TOTAL + DENSTRG (K) *INELASTICPIPA (ATRG (K) ,EPI) 

* 

* end sum over target atoms 

* 

10 CONTINUE 

* 

* 

ELSE IF (J .EQ. 4) THEN 

* 

* sum micro-cross section over target atoms 

* 

DO 20 K=1,NAT 
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* 

* calculation based on scaling from pi p collisions to pi A 

* 

T0TAL= TOTAL + DENSTRG(K) *INELASTICPIMA(ATRG(K) ,EPI) 

* 

* end sum over target atoms 

* 

20 CONTINUE 

* 

* set interaction cross section to zero for muons 

* 

ELSE 

T0TAL=0 . 0 

* 
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END IF 




* add in losses due to decay 




MMLOSS=PIDECAY (EPI , J)+T0TAL 




* end 




RETURN 




* end function mmloss 




SUBROUTINE MULIMITS (EMU , THRESHOLD, UPPERLIM) 


* this subroutine calculates integration limits for muon production 


* from pion decay 
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* 

IMPLICIT NONE 

* 

* Real Variables 

* emu is the muon energy 

* emumev is the muon energy in units of MeV 

* pmu is the muon momentum in umits of Mev 

* threshold is the minimum energy a pion needs to produce a muon 

* of energy emu 

* upperlim is the maximum energy a pion can have and still 

* produce a muon of energy emu 

* pstar moment™ of a produced muon in the pion rest frame 

* estar total energy of a produced muon in the pion 

* rest frame 

* mmu muon mass 

* amu muon mass in atomic mass units 

* mpi pion mass 

* api pion mass in atomic mass units 
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REAL EMU , THRESHOLD , UPPERLIM , EMUMEV , PMUMEV 


$ PSTAR , ESTAR , MMU , MPI , API , AMU 

* 

* START SUBROUTINE MULIMITS 

* 

* set masses 

* 

API= 0.14983 
AMU= 0.11342 
MPI= 139.570 
MMU= 105.658 

* 

* tranform muon energy to units of MeV 

* 

EMUMEV =EMU * AMU 

* 
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* calculate limits 

* 

ESTAR= (MPI**2 . 0+MMU**2 . 0) / (2 . 0*MPI) 

PSTAR=SQRT (ESTAR**2 . 0-MMU**2 . 0) 

PMUMEV=SQRT ( (EMUMEV+MMU) **2 . 0 - MMU**2.0) 

THRESH0LD= (ESTAR* (EMUMEV+MMU) -PSTAR*PMUMEV) *MPI/ (MMU**2 . 0) 

$ - MPI 

UPPERLIM= (ESTAR* (EMUMEV+MMU) +PSTAR*PMUMEV) *MPI/ (MMU**2 . 0) 

$ - MPI 

* 

* convert limits to MeV/amu and prevent mistakes due to round off 

* 

* IF (THRESHOLD .LT. 0.0) THRESHOLDS . 0 
THRESH0LD=THRESH0LD/API 
UPPERLIM=UPPERLIM/ API 

* 
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* end subroutine mulimits 

RETURN 

END 

REAL FUNCTION MUPRODPI (EPI , EMU) 

* This function returns macroscopic cross section for prouction 

* muons from decay of like charged pions 

* 

IMPLICIT NONE 

* 

* Real variables 

* EPI energy of projectile pion that is decaying in MeV/amu 

* EMU energy of muon being created in MeV/amu 

* API mass of pion in MeV/amu 

* AMU mass of muon in MeV/amu 

* ESTAR total energy of muon in pion rest frame 

* PSTAR momentum of muon in pion rest frame 
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* GAMMA lorentz factor of pion relative to medium 

* BETA speed/c of pion relative to medium 

* UPPER upper limit to energy of muon produced from pion of 

* energy EPROJ 

* LOWER lower limit 

* PIDECAY function that returns inverse decay length of pion 

REAL EPI, EMU,EPIMEV, EMUMEV, PSTAR, GAMMA, BETA, ESTAR, 

$ UPPER, LOWER, PIDECAY, API, AMU, MPI, MMU 

* 

* start function muprodpi 

* 

* set masses 

* 

API= 0.14983 
AMU= 0.11342 
MPI= 139.570 
MMU= 105.658 


114 



* 

* tranform muon and pion energy to units of MeV 

* 

EMUMEV=EMU*AMU 

EPIMEV=EPI*API 

* 

* calculate kinematic quantities 

* 

GAMMA= (EPIMEV + MPI)/MPI 
BETA = (1.0 - 1 . 0/ (GAMMA**2 . 0) ) **0 . 5 
ESTAR= (MPI**2 . 0+MMU**2 . 0) / (2 . 0*MPI) 
PSTAR=SQRT (ESTAR**2 . 0-MMU**2 . 0) 
UPPER=GAMMA* (ESTAR+BETA*PSTAR) -MMU 
L0WER=GAMMA* (ESTAR-BETA*PSTAR) -MMU 

* 

* calculate cross section 
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* convert pstar to Mev/amu 

PSTAR=PSTAR/ AMU 

* 

IF ( (EMUMEV .LE. UPPER) .AMD. (EMUMEV .GE. LOWER)) THEM 
MUPRODPI=PIDECAY (EPI , 3) / (2 . 0*BETA*GAMMA*PSTAR) 

ELSE 

MUPR0DPI=0 . 0 
END IF 

* 

* end function muprodpi 

RETURN 

END 

subroutine output (psi ,psimm, st , et) 
implicit none 
integer II ,IJ,IM,ie, j 
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parameter (11=45 , IJ=59, IM=4) 

real psi(ii,ij), psimm(ii , im) , st(ii), et(ii), 

$ tempi, anumm, anu 

open(12,f ile^mmflux.dat’ , status= ; old’ , position^ append’) 
open ( 13 ,f ile= ’ per cent _f lux . dat ’ , status= ’old’ ,position=’ append’ ) 

DO 70 IE=1,II 

write (12 , 120) et(IE), psimm(ie J l)/(anumm(l)*st(ie)) , 

$ psimm(ie , 2) / (anumm(2) *st (ie) ) , 

$ psimm(ie,3)/(anumm(3)*st(ie)) , 

$ psimm(ie,4)/(anumm(4)*st(ie)) , 

$ psimm(ie , 1) / (anumm(l) *st (ie) )+ 

$ psimm(ie , 2) / (anumm(2) *st (ie) )+ 

$ psimm(ie,3)/(anumm(3)*st(ie))+ 

$ psimm(ie,4)/(anumm(4)*st(ie)) 

120 Format (Ell .4, Ell. 4, Ell. 4, Ell. 4, Ell. 4, Ell. 4) 

templ=0 . 0 
do 111 j=l ,2 

templ= templ+psi (ie , J) 
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Ill 


continue 


do 80 J=3, IJ 

templ= templ+psi (ie , J)/anu(J) 

80 continue 

do 90 J=1,IM 

templ=templ+psimm(ie , J)/anumm( J) 

90 continue 

write (13 , 130) et(IE), 100 . *psimm(ie , 1) / (anumm(l) *templ) , 
$ 100.* psimm(ie , 2)/(anumm(2)*templ) , 

$ 100.* psimm(ie ,3)/(anumm(3)*templ) , 

$ 100.* psimm(ie J 4)/(anumm(4)*templ) , 

$ 100.* psimm(ie ,l)/(anumm(l)*templ) + 

$ 100.* psimm(ie ,2)/(anumm(2)*templ)+ 

$ 100.* psimm(ie ,3)/(anumm(3)*templ)+ 

$ 100.* psimm(ie,4)/(anumm(4)*templ) 

130 Format (Ell A, Ell. 4, Ell. 4, Ell. 4, Ell. 4, Ell. 4) 

70 continue 
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close (12) 


close (13) 

return 

end 

* end subroutine output 


SUBROUTINE PGTMM(II, IJ, E, R, ST, PSI, PSIMM, 

$ H, QMOLD) 

* This subprogram propagates the meson and muon components of 

* radiation a distance H. 

* 

IMPLICIT NONE 

* 

* Integer Variables 

* II is the number of range and energy bins 

* IJ is the number of ions 
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* IM is the number of mesons and muons being transported 

* 1 is mu+ , 2 mu-, 3 pi+, 4 pi- 

* I,J are loop indicies 

* IE is loop index for energy and range bins 

INTEGER II, IJ, IM, I, J, IE 
PARAMETER (IM=4) 

* Real Variables 

* r is the range array for protons 

* e is energy in MeV/Amu 

* flux are the fluxes of ions 

* psi are the fluxes of ions multiplied by stopping power 

* fluxmm are the fluxes of mesons and muons 

* psimm fluxes of mesons and muons multiplied by stopping power 

* mmsource is the meson/muon source term 

* qm array of mmsources mutiplied by meson stopping power 

* qmold is qm from the previous iteration 

* deltaqm is the spatial derivitive of qm 

* st is proton stopping power 

* scale stopping power scaling factor needed for source terms 
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* aloss total interaction cross section + inverse decay length 

* qm is the pion and meson source terms 

* H is the stepsize 

REAL E(II) , R(II), PSICII, IJ), FLUX(II , I J) , PSIMM(II , IM) , 

$ FLUXMM(II , IM) , QMOLDdI, IM) , QM(II,IM), ST(II), 

$ MMSOURCE, DELTAQM, SCALE, H, TEMPI, 

$ ALOSS, TEMP 2 , THRESHOLD, UPPERLIM 

* Functions Called 

* anumm calculates stopping power scaling factors for mesons 

* anu calculates stopping power scaling factors for ions 

* mmloss calculates total interaction cross section + decay length 

* phi interpolates flux arrays; in original hzetrn 

* source calculates meson sources at a specific energy 

* mgaussmm integrates source over energy 

* mulimts subroutine caculates limits of integration 

* for muon production 

REAL ANUMM, ANU, MMLOSS, PHI, SOURCE, MGAUSSMM 
EXTERNAL SOURCE 
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* START PROPAGATION 

* 

* calculate fluxes for source term by dividing by stopping power 

* 

DO 10 IE=1,II 
DO 20 J=1,IM 

FLUXMM ( IE , J) =PSIMM ( IE , J) / (ANUMM ( J) *ST ( IE) ) 

20 CONTINUE 

DO 140 J=1 , 2 

FLUX (IE, J)=PSI (IE, J) /ST(IE) 

140 CONTINUE 

DO 100 J=3 , 1 J 

FLUX (IE, J)=PSI(IE, J)/(ANU(J)*ST(IE)) 

100 CONTINUE 

10 CONTINUE 

* 
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* loop for meson and muon particle types 


* 

DO 30 J=1,IM 

* 


* calculate stopping power scaling factor 


* 

SCALE=ANUMM ( J ) 

* 


* loop for range grid 


* 

DO 40 IE=1,II 

* 


* calculate rate of particle loss at energy IE 


* 
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ALOSS= MMLOSS (E (IE) , J) 


* 

* calculate integration limits for particle production 

* Keep them within the energy range being transperted 

* 

IF ((J.EQ.l) .OR. (J.EQ.2)) THEM 

CALL MULIMITS(EQE) , THRESHOLD , UPPERLIM) 
IF (THRESHOLD .LT. E(l)) THRESHOLD=E ( 1) 
IF (UPPERLIM .GT. E(II)) UPPERLIM=E (II) 
ELSE 

* IF (E (IE) .GT. 290.0) THEN 

* THRESHOLD=E (IE) 

* ELSE 

THRESH0LD=290 . 0 

* END IF 
UPPERLIM=E(II) 

END IF 

* print*, threshold, upperlim, E(IE) 

* 
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* calculate source terms 

* 

MMSOURCE=MGAUSSMM (THRESHOLD , UPPERLIM , 6 , SOURCE , J , 

$ E(IE) , E , FLUX , FLUXMM , II , I J , IM) 

QM(IE, J)=SCALE * ST (IE) *MMSOURCE 
DELTAQM= ( QM ( IE , J) -QMOLD ( IE , J ) ) /H 

* 

* calculate new psimm 

* note that (h/2) * (1 . 0+0 . 66666666*a*h +((a*h)**2.0)*exp(ah)/12.0 

* = exp(ah)*(h/a-l/a**2) +l/a**2 

* 

TEMP1=QM(IE , J) * (EXP (AL0SS*H) -1 . 0) /ALOSS 
TEMP 2= (H**2 . 0/2 . 0) * (1 . 0+0 . 66666666*AL0SS*H 
$ + ( (AL0SS*H) **2 . 0) *exp(AL0SS*H)/12 . 0) 

PSIMM (IE, J)=EXP (-H+AL0SS) * 

$ (PHI(R(II)+1.0,II,R,PSIMM(1, J) ,R(IE) + SCALE*H) 

$ +TEMP1 + DELTAQM*TEMP2) 
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* 

* end loop over range 

* 

40 CONTINUE 

* 

* end loop over meson type 

* 

30 CONTINUE 

* 

* set source term for use in next iteration 

* 

DO 50 J=1,IM 
DO 60 IE=1 , II 

QM0LD(IE, J)=QM(IE, J) 

60 CONTINUE 
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50 CONTINUE 


* 

* end subroutine pgtmm 

RETURN 

END 

REAL FUNCTION PIDECAY(EPI, J) 

* This function returns the inverse decay length of pions/muons of 

* type J and energy EPI. 

* 

IMPLICIT NONE 

* 

* Integer Variables 

* IM number of mesons and muons 

* J indexes particle type 1, 2, 3 or 4 correspond to 
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mu+, mu-, pi+, and pi- 


* 

INTEGER IM, J 
PARAMETER (IM=4) 

* Real Variables 

* EPI is the energy in MeV/amu 

* EMEV is the energy in MeV 

* TIME is the lifetime in rest frame in seconds 

* C is the speed of light in cm/s 

* MPI are particle masses in MeV/c~2 

* API is mass in amu 

* GAMMA is the time dilation factor. 

* BETA is the velocity/c 

* DN is the target density in gm/cnT3 

REAL EPI, TIME(IM) , C, MPI(IM), GAMMA, BETA, API(IM), EMEV, DN 
PARAMETER (C=2 . 9979245800E10) 

DATA TIME /2.19703E-6, 2.19703E-6, 2.6033E-8, 2.6033E-8/ 

DATA MPI /105.658, 105.658, 139.57, 139.57/ 

DATA API /0. 11342, 0.11342, 0.14983, 0.14983/ 
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* unused reals for common statement 

REAL NLAY , XLAY , MAT , ZTRG ( 5 ) , ATRG ( 5 ) , DENSTRG ( 5 ) 

* set commons 

COMMOM/TARGET/MLAY , XLAY , DM , MAT , ATRG , ZTRG , DEMSTRG 

* trans form energy from MeV/amu to MeV 

* 

EMEV=EPI*API ( J) 

* 

* calculate inverse decay length 

* 

GAMMA= (EMEV + MPI (J))/MPI (J) 

BETA = (1.0 - 1 . 0/ (GAMMA**2 . 0) ) **0 . 5 
PIDECAY=1 . 0/ (BETA*C*TIME(J) *GAMMA*DN) 
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* 

* end function pidecay 

RETURN 

END 

REAL FUNCTION PIMAPIM(EPROJ ,EPI , AT) 

* this function that returns microscopic spectral distribution for 

* the reaction pi- + A -> pi- + x 

* 

IMPLICIT NONE 

* 

* Real Variables 

* EPROJ kinetic energy of projectile pi+ 

* EPI kinetic energy of produced pi+ 

* PIMINUS function returns spectral distribution for reaction 
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proton + proton -> pi- +x 


* 

* ALPHA scales collisions from target nucleon to nucleus 

* AT atomic number of target 

REAL EPRO J , EPI , PIMINUS, ALPHA , AT 

* 

start function PIMAPIM 

* set alpha 

* 

ALPHA =0.77 

* 

* calculate cross section 

* this formulae assumes pion projectiles act like protons 

* 

PIMAPIM=0 . 62*PIMINUS (EPRO J , EPI) *AT**ALPHA 
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* 

* end function pimapim 

RETURN 

END 

REAL FUNCTION PIMAPIP (EPROJ ,EPI , AT) 

* this function that returns microscopic spectral distribution for 

* the reaction pi- + A -> pi+ + x 

* 

IMPLICIT NONE 

* 

* Real Variables 

* EPROJ kinetic energy of projectile pi- 

* EPI kinetic energy of produced pi+ 

* PIPLUS function returns spectral distribution for reaction 
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proton + proton -> pi+ +x 


* 

* ALPHA scales collisions from target nucleon to nucleus 

* A atomic number of target 

REAL EPRO J , EPI, PIPLUS, ALPHA , AT 

* 

start function PIMAPIP 

* set alpha 

* 

ALPHA =0.77 

* 

* calculate cross section 

* this formulae assumes pion projectiles act like protons 

* 

PIMAPIP=0 . 62*PIPLUS (EPRO J , EPI) *AT**ALPHA 
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* 

* end function pimapip 

RETURN 

END 

REAL FUNCTION PIMINUS(EPR, E) 

* This function returns microscopic spectral distribution for 

* negative pion production from proton-proton collisions using a 

* parameterization. The formula is split into three different 

* energy regions 

* 

IMPLICIT NONE 

* 
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* Real Variables 


* G and H parameters of parametrization 

* EPR energy of projectile proton in MeV 

* E energy of pion in MeV 

* EPROJ energy of projectile proton in GeV 

* EPION energy of pion in GeV 

REAL TEMPI, TEMP 2, G(14), H(17), EPROJ, EPI, EPR, E, API 

start function PIMINUS 

* set parameters 

* 

DATA G / 1 . 06E-9, -2 .8 ,3 .7E-8, -1 .89,22.3,-1 .5,-30 .5,0. 0938,-24 .6, 
$ 0. 0313, 0.25, 2. 5E6, 7. 96, -49. 5/ 

DATA H /2 . 39E-10 ,-2.8,1. 14E-8 ,-2.3,22.3,-2.23, -31 .3,0. 0938 , 

$ -24. 9, 0.0313, 2. 5E6, 0.25, 60322, 1.1, -65. 9, -9. 39, -1.25/ 

API=0. 14983 
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* 

* change energy from MeV/amu to GeV, because the parametrization 

* assumes energy is in units of GeV. A proton has an atomic number 

* of one so the amu part does not need to be converted. 

* 

EPI= E*API/1000 . 0 
EPRO J=EPR/ 1000 . 0 

* 

* if pion energy too low set it to prevent overflow 

* 

IF (EPI .LT. 0.000001) THEM 
EPI=0. 000001 
END IF 

* 

* calculate parameterization which is broken up into 3 energy ranges. 
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IF (EPROJ .GE. .29 .AND. EPROJ .LE. 2.0) THEN 


* 


TEMP 2 = G ( 1) *EPI**G(2) + G (3) *EPRO J**G (4) 

TEMP1=EXP (G (5) + G (6) /SQRT (EPROJ) + G(7) *EPI**G(8) + 
$ G(9)*EPI**G(10)) 

PIMINUS = EPI**G(11)*(G(12)*TEMP1/TEMP2 + 

$ G(13)*EXP( G(14)*SQRT(EPI))) 


ELSE IF (EPROJ .GT. 2.0) THEN 


TEMP 2 = H(l) *EPI**H(2) +H (3) *EPRO J**H (4) 

TEMPI = EXP(H(5) + H (6) /SQRT (EPROJ) + H(7) *EPI**H(8) 

$ +H (9) *EPI**H (10) ) 

PIMINUS = H(ll) *EPI**H (12) *TEMP1/TEMP2 + H(13) 

$ *EPI**H(14)*EXP(H(15)*SQRT(EPI)+H(16)*EPR0J**H(17)) 


ELSE 

PIMINUS =0.0 
RETURN 
END IF 
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* 

* change spectral distribution from mb/GeV to cnr2*amu/MeV 

* 

PIMINUS=PIMINUS*API*1 . 0E-30 

* 

* end function piminus 

RETURN 

END 

REAL FUNCTION PIPAPIM(EPROJ,EPI, AT) 

* this function that returns microsopic spectral distributiuon for 

* the reaction pi+ + A -> pi- + x 

* 
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IMPLICIT NONE 


* 

* Real Variables 

* EPROJ kinetic energy of projectile pi+ 

* EPI kinetic energy of produced pi+ 

* PIMINUS function returns spectral distribution for reaction 

* proton + proton -> pi- +x 

* ALPHA scales collsions from target nucleon to nucleus 

* AT atomic number of target 

REAL EPROJ , EPI , PIMINUS, ALPHA, AT 

* 

* start function PIPAPIM 

* 

* set alpha 

* 

ALPHA =0.77 
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* 

* calculate cross section 

* this formulae assumes pion projectiles act like protons 

* 

PIPAPIM=0 . 62*PIMINUS (EPRO J , EPI) *AT**ALPHA 

* 

* end function pipapim 

RETURN 

END 


REAL FUNCTION PIPAPIP (EPRO J , EPI , AT) 

* this function that returns microsopic spectral distributiuon for 

* the reaction pi+ + A -> pi+ + x 
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IMPLICIT NONE 


* 

* 

* Real Variables 

* EPROJ kinetic energy of projectile pi+ 

* EPI kinetic energy of produced pi+ 

* PIPLUS function returns spectral distribution for reaction 

* proton + proton -> pi+ +x 

* ALPHA scales collsions from target nucleon to nucleus 

* AT atomic number of target 

REAL EPROJ, EPI, PIPLUS, ALPHA, AT 

* 

* start function PIAPIP 

* 

* set alpha 

* 


141 




ALPHA =0.77 


* 

* calculate cross section 

* this formulae assumes pion projectiles act like protons 

* 

PIPAPIP=0 . 62*PIPLUS (EPROJ ,EPI) *AT**ALPHA 

* 

* end function pipapip 

RETURN 

END 

REAL FUNCTION PIPLUS(EPR, EPI) 

* This function returns microscopic spectral distribution for 

* positive pion production from proton-proton collisions using a 

* parameterization. The formula is split into three different 
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* energy regions . 

* 

IMPLICIT NOME 

* 

* Real Variables 

* C and D parameters of parameterization 

* EPR energy of projectile proton in MeV 

* EPI energy of pion in MeV 

* EPROJ energy of projectile proton in GeV 

* EPION energy of pion in GeV 

REAL TEMPI , TEMP 2, C(20), D(17), EPROJ, EPION, EPI, EPR, API 

start function PIPLUS 

* 

* set parameters 
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DATA C /2.2E-8, -2. 7,4. 22E-7.-1. 88,22. 3,1. 98,-0. 28,-1. 75, 


* 


$ -29 .4,0. 0938,-24.4,0. 0312,0. 0389, 1.78,2. 5E6 ,0 . 25 , 

$ 976,2.3, -46 , -0 . 989/ 

DATA D /4 . 5E-11 ,-2 .98,1. 18E-9 ,-2 .55,22.3,-0 .765,-35 .3,0. 0938 , 
$ -22. 5, 0.0313, 2. 56E6, 0.25, 60322, 1.18, -72. 2, 0.9421, 0.1/ 


API=0. 14983 

* 

* change energy from MeV/amu to GeV. the parameterization assumes 

* energy is in units of GeV. Note that a proton has an atomic number 

* of one so the amu part does not need to be converted. 

* 

EPI0N=EPI*API/ 1000 . 0 
EPR0 J=EPR/ 1000 . 0 

* 

* if pion energy too low set it to prevent overflow 
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IF (EPION .LT. O.OOOOOl) THEM 


* 


EPION = 0.000001 
END IF 

* 

* calculate parameterization which is broken up into 3 energy ranges. 

* 

IF ((EPROJ .GE. .29) .AMD. (EPROJ .LE. 2.0)) THEM 
TEMP 2 = C (1) *EPI0N**C (2) + C (3) *EPR0 J**C (4) 

TEMP1=EXP (C (5) +C (6) /SQRT (EPROJ) +C (7) *EPR0 J**C (8) 

$ +C(9)*EPI0N**C(10)+C(11)*EPI0N**C(12)*EPR0J**C(13) 

$ +C(14) *L0G (EPROJ) ) 

PIPLUS = C(15)*EPI0N**C(16)*TEMP1/TEMP2 + C(17)* 

$ EPI0N**C(18)*EXP(C(19)*SQRT(EPI0N) 

$ + C (20) *SQRT (EPROJ) ) 

ELSE IF (EPROJ .GT. 2.0) THEM 

TEMP 2 = D (1) *EPION**D (2) +D (3) *EPRO J**D (4) 
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TEMPI = EXP (D (5) + D (6) / SQRT (EPRO J) + D (7) *EPION**D (8) 
$ +D (9) *EPION**D (10) ) 

PIPLUS = D(11)*EPI0N**D(12)*TEMP1/TEMP2 + 

$ D ( 13) *EPION**D (14) *EXP (D (15) *SQRT (EPION) 

$ +D (16) *EPROJ**D (17) ) 

ELSE 

PIPLUS =0.0 
END IF 

* 

* change spectral distribution from mb/GeV to cnT2*amu/MeV 

* lmb=10~-27 cnT2 

* 

PIPLUS=PIPLUS*1 . 0E-30*API 

* 

* end function piplus 
RETURN 
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END 


REAL FUNCTION SOURCE (EPI, EPROJ, E, FLUX, FLUXMM , J, II, 
$ IJ, IM) 

* This function calculates meson sources at specific energies. 

* It is integrated over that energy by mgaussmm. After this 

* function is integrated by mgaussmm, it is the function 

* q_j(x,r_i) described in the text. 

* 

IMPLICIT NONE 

* 

* Integer Variables 

* II number of range grid points 

* IJ total number of ions 

* IM total number of mesons/muons 

* P index for projectiles 

* J type of meson being produced 
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INTEGER II, IJ, P, J, IM 

* Real Variables 

* E energy grid 

* FLUX ion flux array 

* FLUXMM meson flux array 

* EPI energy of meson being produced 

* EPROJ energy of projectile 

* FLUENCE interpolated value of flux 

* AAMMCS function that returns macroscopic cross section for 

* production of mesons from collisions of nuclei with 

* the medium 

* MAMMCS function that returns macroscopic cross section for 

* production of mesons from collisions of mesons with 

* the medium 

* MUPRODPI function that returns macroscopic cross section for 

* production muons from decay of like charged pions 

* PHI function that interpolates fluxes 

* TEMPSOURCE 
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REAL E(II), FLUX(II , I J) , FLUXMM(II , IM) , EPI , EPRO J , FLUENCE 


$ AAMMCS, MAMMCS , PHI , MUPRODPI, TEMPSOURCE 

* 

* start function 

* 

* initialize source 

* 

TEMPSOURCE=0 . 0 

* 

* Ignore muon production from collisions 

* 

IF (J .GT. 2) THEM 

* 

* sum over projectile nucleons 
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DO 50 P=1,IJ 


* 

* 

* calculate projectile flux at energy eproj . 

* 

FLUENCE=PHI (E (II) +1 . 0 , II , E , FLUX ( 1 , P) , EPROJ) 

TEMPS 0URCE=TEMP SOURCE+FLUENCE* AAMMC S (EPROJ , EPI , P , J) 

* 

* end sum over projectile nucleons 

* 

50 CONTINUE 

* 

* sum over meson projectiles, no muons included 

* 

DO 60 P=3 , IM 
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* 

* calculate projectile flux at energy eproj . 

* 

FLUENCE=PHI (E(II)+1.0,II,E, FLUXMM ( 1 , P ) , EPRO J ) 
TEMPSOURCE=TEMPSOURCE+FLUENCE*MAMMCS (EPROJ , EPI , P , J) 

* 

* end sum over meson projectiles 

* 

60 CONTINUE 

* 

* now do muon production 

* muons are being produced only from pion decay. 

* kaon are not yet included. 

* 

ELSE 
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FLUENCE=PHI (E(II)+1.0,II,E, FLUXMM ( 1 , J+ 2 ) , EPRO J ) 


TEMPSOURCE=TEMPSOURCE+FLUENCE*MUPRODPI (EPROJ, EPI) 

* 

* end function source 
END IF 

SOURCE=TEMPSOURCE 

RETURN 

END 
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